source/ga/GA.f

   1       program unresag
   2
   3       implicit none
   4       include 'GA.inc'
   5       include 'common.inc'
   6       include 'io.inc'
   7       real*8 :: WEIGHTS(18)
   8       integer*4 :: i,j,m,n,iloc,pos,iter,maxiter,numarg,tmplen
   9       integer*4 :: WybierzOsobnika,Najlepszy,ZnajdzPodobnego,FindWorst
  10       integer*4 :: status
  11       integer*4 :: mnoznik = 10
  12       character*16 :: mode,argn,argi,dir_nam
  13       character*15 :: filename
  14       real*8 :: loc,fc,sumfitness,sump,avrp,maxfc
  15       real :: rand, r
  16       real*8 :: maxfcever
  17       real*8 :: FunkcjaCeluB, FunkcjaCeluB2, FunkcjaCeluR
  18       character(80),dimension(5) :: wagi
  19       logical :: searching = .false.
  20       logical :: logdata = .true.
  21       logical :: fileex = .false.
  22       logical :: retry = .false.
  23       logical :: dir_ex = .false.
  24       logical :: debug = .true.
  25       character*200 :: text, tmptext,tmptext2
  26       character(len=*), parameter :: FMT = "(19F10.5,E15.7,F10.5)"
  27
  28 c =======================================================================
  29 c  Main program
  30 c =======================================================================
  31
  32
  33       call ReadInput(status)
  34       if ((maxminstep.gt.0).and.(generation.eq.0)) then
  35 c Do first score
  36        do_fs=.true.
  37       endif
  38 c  Is everything  OK?
  39       if (status.gt.0) then
  40        write(*,*) "There ware input file errors - check log."
  41        goto 2010
  42       endif
  43
  44       call itime(timeArray)
  45       r = rand(timeArray(1)*10000+timeArray(2)*100+timeArray(3))
  46 c     Fair enough?
  47 c      r = 527865426
  48
  49
  50 c-------------------------------------------------------
  51 c bank & populacja arrays store this
  52 c 1-19 - weights
  53 c   20 - energy obtained from zscore
  54 c   21 - fittness
  55 c-------------------------------------------------------
  56       allocate (bank(banksize,21))
  57       allocate (populacja(BANK_MULTIPLIER*banksize,21))
  58       allocate (pairs(BANK_MULTIPLIER*banksize))
  59 c
  60 c Zero the arrays
  61 c
  62       do i=1,BANK_MULTIPLIER*banksize
  63        do j=1,21
  64         populacja(i,j)=0
  65        enddo
  66        pairs(i)=0
  67       enddo
  68       do i=1,banksize
  69        do j=1,21
  70        bank(i,j)=0
  71        enddo
  72       enddo
  73
  74 c First time here?
  75       inquire(FILE=ostatefn,EXIST=fileex)
  76 c No. Prepere next generation
  77       if (fileex) then
  78        open(ostate,file=ostatefn)
  79        read(ostate,'(I4)') generation
  80        read(ostate,'(L2)') do_optima
  81        read(ostate,'(L2)') do_ga
  82        read(ostate,'(F7.5)') avrd
  83        close(ostate)
  84        if (do_ga) then
  85         call write2log("Doing GA in this step")
  86        endif
  87        if (.not.do_optima) then
  88         maxminstep=0
  89         call write2log("ZSCORE weights minimalization disabled for now")
  90         generation=generation+1
  91        endif
  92        write(tmptext,'(I4)') generation
  93        call write2log("This is genaration "//tmptext)
  94        call ReadOptimaW(BANK_MULTIPLIER*banksize,populacja)
  95
  96 c Yes. Generate random zero-population
  97       else
  98        call GenPopulation(BANK_MULTIPLIER*banksize,populacja)
  99        write(tmptext,'(I4)') generation
 100        call write2log("This is genaration "//tmptext)
 101        if (do_optima) then
 102         do_ga=.false.
 103        endif
 104
 105 c End of "fisrst time here?" code
 106       endif
 107
 108 c
 109 c Do we actual use a Genetic Algorithm?
 110 c
 111
 112       if (do_ga.eqv..true.) then
 113
 114 c Yes. This is only done when doing second pass with ZSCORE (without weight minimalizaton)
 115 c or with weight minimalization disabled (maxmin=0) just to obtain final score)
 116
 117 c \\---//================================================================c
 118 c  \\-//                                                                 c
 119 c   "//  Genetic Algorithm code starts here                              c
 120 c   //-\\                                                                c
 121 c =//---\\===============================================================c
 122
 123       call ReadZEnergy(BANK_MULTIPLIER*banksize,populacja)
 124       do i=1,BANK_MULTIPLIER*banksize
 125 c       populacja(i,20)=FunkcjaCeluR(populacja(i,20))
 126        write(*,*) "Ind ",i,populacja(i,20)
 127       enddo
 128
 129       call write2log("==[ Doing GA now ]=============================")
 130
 131 c
 132 c Populate bank depending on algorithm
 133 c
 134       select case(alg)
 135        case('simple')
 136         call GetNBest(populacja,bank,banksize)
 137        case('better')
 138         call GetNBest(populacja,bank,banksize)
 139        case('betterv2')
 140         call GetNBest(populacja,bank,banksize)
 141        case('csa')
 142 c --- debug begin ---
 143         if (debug) then
 144          write(*,*) "generation: ",generation
 145          write(*,*) "maxminstep: ",maxminstep
 146          write(*,*) "do_ga: ",do_ga
 147          write(*,*) "do_optima: ",do_optima
 148          write(*,*) "do_fs: ",do_fs
 149 c        write(*,*) "cicutoff: ", cicutoff
 150         endif
 151 c --- debug end ---
 152
 153 c        csacutoff=cicutoff
 154
 155 c
 156 c  Fill the bank just after the first time we get the score
 157 c
 158         if (((generation.eq.1).and.(maxminstep.eq.0)).or.((generation.eq&
 159      &.1).and.(maxminstep.gt.0))) then
 160
 161 c --- debug begin ---
 162          if (debug) then
 163           write(*,*) "Calling GetNBest..."
 164          endif
 165 c --- debug end ---
 166          call GetNBest(populacja,bank,banksize)
 167          call CalcAvgDist(bank,avrd)
 168          write(tmptext,'(F7.5)') avrd
 169          call write2log("Average distance between individuals in initial&
 170      & bank is "//trim(tmptext))
 171          csacutoff=(maxco*avrd)-generation*avrd*(maxco-minco)/maxgen
 172          write(tmptext,'(F7.5)') csacutoff
 173          call write2log("CSA cutoff is now set to "//trim(tmptext))
 174 c         csacutoff=maxco*avrd
 175
 176
 177
 178 c --- debug begin ---
 179          if (debug) then
 180           do i=1,banksize
 181            write(*,FMT) bank(i,:)
 182           enddo
 183          endif
 184 c --- debug end ---
 185
 186          call CalcFitness(bank,sumfitness)
 187          call WriteBank(bank)
 188
 189 c
 190 c Do normal CSA
 191 c
 192         else
 193
 194 c --- debug begin ---
 195          if (debug) then
 196           write(*,*) "Szukam podobnego"
 197          endif
 198 c --- debug end ---
 199
 200          call ReadBank(bank)
 201          write(tmptext,'(F7.5)') avrd
 202          call write2log("Average distance in bank is "//trim(tmptext))
 203
 204          csacutoff=maxco*avrd-generation*avrd*(maxco-minco)/maxgen
 205          write(tmptext,'(F7.5)') csacutoff
 206          call write2log("CSA cutoff is now set to "//trim(tmptext))
 207
 208          do i=1,BANK_MULTIPLIER*banksize
 209           write(tmptext,'(I4)') i
 210           call write2log("Checking ind "//trim(tmptext))
 211           j=ZnajdzPodobnego(banksize,bank,populacja(i,:),csacutoff)
 212           if (j.gt.0) then  ! W banku jest podobny
 213            if (populacja(i,20).lt.bank(j,20)) then
 214             write(tmptext,'(I4)') j
 215             write(tmptext2,'(I4)') i
 216             call write2log("  Swaping ind"//trim(tmptext)//" from bank t&
 217      &o ind "//trim(tmptext2)//" from population")
 218             write(tmptext2,'(19F8.5,E15.7,F8.5)') bank(j,:)
 219             call write2log("  BANK"//trim(tmptext)//":"//trim(tmptext2))
 220
 221             write(tmptext,'(I4)') i
 222             write(tmptext2,'(19F8.5,E15.7,F8.5)') populacja(i,:)
 223             call write2log("  POP "//trim(tmptext)//":"//trim(tmptext2))
 224
 225             bank(j,:)=populacja(i,:)
 226            else
 227             call write2log("  Found simialar but not better")
 228            endif
 229           else              ! W banku nie ma podobnego
 230            j=FindWorst(banksize,bank)
 231            write(tmptext2,'(I4)') j
 232            if (populacja(i,20).lt.bank(j,20)) then
 233             call write2log("  Worst in bank is "//trim(tmptext2))
 234             write(tmptext,'(I4)') i
 235             call write2log("  Swaping worst ind in bank to "//trim(tmpte&
 236      &xt))
 237             write(tmptext,'(I4)') j
 238             write(tmptext2,'(19F8.5,E15.7,F8.5)') bank(j,:)
 239             call write2log("  BANK"//trim(tmptext)//":"//trim(tmptext2))
 240             write(tmptext,'(I4)') i
 241             write(tmptext2,'(19F8.5,E15.7,F8.5)') populacja(i,:)
 242             call write2log("  POP "//trim(tmptext)//":"//trim(tmptext2))
 243             bank(j,:)=populacja(i,:)
 244            else
 245             call write2log("  The worst in bank is better then this Ind"&
 246      &)
 247            endif
 248           endif
 249          enddo
 250
 251 c -----------------------------------------------------------------------
 252 c      Calculate fitness for reproduction
 253 c -----------------------------------------------------------------------
 254           call CalcFitness(bank,sumfitness)
 255
 256 c --- debug begin ---
 257         if (debug) then
 258          write (*,*) "Dumping bank: "
 259          do i=1,banksize
 260           write(*,FMT) bank(i,:)
 261          enddo
 262          write (*,*) "Sumfitness: ", sumfitness
 263         endif
 264 c --- debug end ---
 265
 266          call WriteBank(bank)
 267
 268         endif
 269        case('cluster')
 270         write(*,*) "Some stuff here in the future"
 271         goto 2010
 272       end select
 273
 274 c      call WritePopSum()
 275
 276
 277 c Have we done last generation ?
 278       if (generation.eq.maxgen) then
 279        call write2log("Maximum number of genarations reached. QUITING.")
 280        goto 2010
 281       endif
 282 c Prapering for  next generation so increase counter
 283         iter=iter+1
 284
 285 c -----------------------------------------------------------------------
 286 c       Reproduction
 287 c -----------------------------------------------------------------------
 288       call write2log("==[ Reproducting ]============================")
 289       select case(alg)
 290        case('simple','csa')
 291         do i=1,BANK_MULTIPLIER*banksize
 292         m=WybierzOsobnika(banksize,bank)
 293         write(tmptext,'(I4)') m
 294         write(tmptext2,'(I4)') i
 295         call write2log("Reproducting individual "//trim(tmptext)//" from&
 296      & bank as new individual "//trim(tmptext2))
 297         populacja(i,:)=bank(m,:)
 298         end do
 299        case('better','betterv2')
 300         populacja(1,:)=bank(1,:)
 301         call write2log("Reproducting first (best) individual")
 302         do i=2,BANK_MULTIPLIER*banksize
 303         m=WybierzOsobnika(banksize,bank)
 304         write(tmptext,'(I4)') m
 305         write(tmptext2,'(I4)') i
 306         call write2log("Reproducting individual "//trim(tmptext)//" from&
 307      & bank as new individual "//trim(tmptext2))
 308         populacja(i,:)=bank(m,:)
 309         end do
 310       end select
 311
 312 c -----------------------------------------------------------------------
 313 c       Pairing
 314 c -----------------------------------------------------------------------
 315       do i=1,BANK_MULTIPLIER*banksize
 316        pairs(i)=0
 317       end do
 318
 319       do i=1,BANK_MULTIPLIER*banksize
 320        if (pairs(i).eq.0) then
 321 210     pos=1+int(rand(0)/(1/real(BANK_MULTIPLIER*banksize)))
 322         if (pos.le.i) then
 323          goto 210
 324         endif
 325         do j=1,(i-1)
 326          if (pairs(j).eq.pos) then
 327           goto 210
 328          endif
 329         end do
 330         pairs(i) = pos
 331         pairs(pos) = i
 332        endif
 333       end do
 334       call write2log("==[ Pairing ]==============================")
 335       do i=1,BANK_MULTIPLIER*banksize
 336        write(tmptext,'(A11,I3.3,A5,I3.3,A21)') "Individual ",i," and ",p&
 337      &airs(i)," will be crossed over"
 338        call write2log(trim(tmptext))
 339       end do
 340
 341 c -----------------------------------------------------------------------
 342 c       Crossover
 343 c -----------------------------------------------------------------------
 344       call write2log("==[ Crossing ]==============================")
 345       do i=1,BANK_MULTIPLIER*banksize
 346        if (pairs(i) .gt. i) then
 347         call CrossoverS(populacja(i,:),populacja(pairs(i),:))
 348        endif
 349       end do
 350
 351
 352        call write2log("         WLONG   WSCP    WELEC   WBOND   WANG    &
 353      &WSCLOC  WTOR    WTORD   WCORRH  WCORR4  WCORR5  WCORR6  WEL_LOC WT&
 354      &URN3  WTURN4  WTURN6  WVDWPP  WHPB    WSCCOR")
 355       do i=1,BANK_MULTIPLIER*banksize
 356        tmplen = 9
 357        text = ""
 358        write(text,'(A4,I3.3,A2)') "Ind",i,": "
 359        do j=1,19
 360         write(tmptext,'(F7.5)') populacja(i,j)
 361         text =  text(1:tmplen) // tmptext
 362         tmplen = tmplen + 8
 363        end do
 364        call write2log(trim(text))
 365       end do
 366
 367 c -----------------------------------------------------------------------
 368 c       Mutation
 369 c -----------------------------------------------------------------------
 370
 371       call write2log("==[ Mutating ]==============================")
 372       do i = 1,19
 373        if (IGNORE_WEIGHT(i).eq.1) then
 374         write(tmptext,'(I2)') i
 375         call write2log("Parameter "//trim(WNAME(i))//" will be ignored")
 376        endif
 377       end do
 378       do i=1,BANK_MULTIPLIER*banksize
 379        call MutationV(populacja(i,:),i,MUTRATIO)
 380       end do
 381 c for betterv2 and csa restore best individual
 382       select case(alg)
 383        case('betterv2','csa')
 384        populacja(1,:)=bank(1,:)
 385       end select
 386
 387 c
 388 c Genetic Algorithm blokcode stops here
 389 c
 390       endif
 391
 392 c ---------------------------------
 393 c  Do some final stuff
 394 c ---------------------------------
 395 c
 396 c Set the variables before writing to state file
 397 c
 398       if ((.not.do_ga).and.(.not.do_optima)) then
 399        if (do_fs) then
 400         do_ga=.false.
 401         do_optima=.true.
 402        else
 403         do_ga=.true.
 404         do_optima=.false.
 405        endif
 406       else
 407        if (do_fs) then
 408         do_optima=.not.do_optima
 409         if (generation.eq.0) then
 410          do_ga=.false.
 411         else
 412          do_ga=.not.do_ga
 413         endif
 414        else
 415         do_ga=.true.
 416         do_optima=.false.
 417        endif
 418       endif
 419
 420
 421 c
 422 c Write te current state and population
 423 c
 424       call WriteState()
 425 c      call WritePopSum()
 426 c
 427 c Create the inputs
 428 c
 429       write(tmptext,'(I3)') generation+1
 430       call write2log("Preparing inputs for next generation ("//trim(adju&
 431      &stl(tmptext))//")")
 432       call CreateInputs(BANK_MULTIPLIER*banksize,populacja)
 433 c
 434 c All done? Then let's quit.
 435 c
 436       goto 2010
 437
 438
 439
 440 c ======================================================================
 441 c
 442 c    MAIN CODE ENDS HERE
 443 c
 444 c ======================================================================
 445
 446
 447
 448 2000  deallocate(populacja)
 449       deallocate(pairs)
 450 2010  end program
 451
 452
 453 c ======================================================================
 454 c  Clustering subroutines
 455 c ======================================================================
 456
 457       include 'cluster.inc'
 458
 459
 460 c ======================================================================
 461 c  Write2log subroutine
 462 c ======================================================================
 463 c  Creates logs
 464 c ----------------------------------------------------------------------
 465
 466       subroutine Write2log(text)
 467        include 'common.inc'
 468        include 'io.inc'
 469        character*200 :: tmptext
 470        character(len=*) :: text
 471        logical :: ex
 472        call itime(timeArray)
 473        inquire(file=ologfn,exist=ex)
 474        if (.not.ex) then
 475        open(olog, file=ologfn, access="append")
 476        write(olog,"(A)") "==============================================&
 477      &=========="
 478        tmptext = "= UNRES weights optimizer version "//version
 479        write(olog,"(A)") tmptext
 480        write(olog,"(A)") info
 481        write(olog,"(A)") "==============================================&
 482      &=========="
 483        close(olog)
 484        endif
 485        open(olog, file = ologfn , access = "append")
 486        write(olog,"(I2.2,A1,I2.2,A1,I2.2,A2,A)") timeArray(1),":",timeAr&
 487      &ray(2),":",timeArray(3),": ",text
 488        close(olog)
 489       end subroutine Write2log
 490
 491 c ======================================================================
 492 c  GetNBest subroutine
 493 c ======================================================================
 494 c Fills up the bank population up to banksize with individuals from
 495 c inputpop having the best(lowest) score
 496 c ----------------------------------------------------------------------
 497
 498       subroutine GetNBest(inputpop,bank,banksize)
 499        integer*4 :: banksize, i,j,idx
 500        real*8 :: inputpop, bank, best, last
 501        include 'GA.inc'
 502        dimension inputpop(banksize*BANK_MULTIPLIER,21)
 503        dimension bank(banksize,21)
 504
 505        last=0.0
 506        idx=1
 507        do j=1,banksize
 508        best=huge(0.0d0)
 509         do i=1,banksize*BANK_MULTIPLIER
 510          if ((inputpop(i,20).lt.best).and.(inputpop(i,20).gt.last)) then
 511           best=inputpop(i,20)
 512           idx=i
 513          endif
 514         end do
 515         bank(j,:)=inputpop(idx,:)
 516         last=best
 517        end do
 518       end subroutine GetNBest
 519
 520 c ======================================================================
 521 c  FindWorst subroutine
 522 c ======================================================================
 523 c  Finds the worst (highest score) individual in population
 524 c ----------------------------------------------------------------------
 525
 526       integer*4 function FindWorst(rozmiar,pop)
 527        integer*4 :: rozmiar, i, idx
 528        real*8 :: pop,last
 529        dimension pop(rozmiar,21)
 530
 531        last=0.0
 532        idx = 0
 533        do i=1,rozmiar
 534         if (pop(i,20).gt.last) then
 535          idx = i
 536          last=pop(i,20)
 537         endif
 538        end do
 539        FindWorst=idx
 540        return
 541       end function
 542
 543 c ======================================================================
 544 c  GenPopulation subroutine
 545 c ======================================================================
 546 c  Generates a population (pop) of nind individuals with 19 random
 547 c  wights from a set of ranges (see GA.inc).
 548 c ----------------------------------------------------------------------
 549
 550       subroutine GenPopulation(nind,pop)
 551        integer*4 :: i,j,nind
 552        real*8 :: pop
 553        dimension pop(nind,21)
 554        include 'GA.inc'
 555        do j=1,nind
 556         do i=1,19
 557          select case(i)
 558           case (1)
 559            pop(j,1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
 560           case (2)
 561            pop(j,2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
 562           case (3)
 563            pop(j,3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
 564           case (4)
 565            pop(j,4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
 566           case (5)
 567            pop(j,5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
 568           case (6)
 569            pop(j,6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
 570           case (7)
 571            pop(j,7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
 572           case (8)
 573            pop(j,8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
 574           case (9)
 575            pop(j,9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
 576           case (10)
 577            pop(j,10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
 578           case (11)
 579            pop(j,11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
 580           case (12)
 581            pop(j,12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
 582           case (13)
 583            pop(j,13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
 584           case (14)
 585            pop(j,14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
 586           case (15)
 587            pop(j,15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
 588           case (16)
 589            pop(j,16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
 590           case (17)
 591            pop(j,17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
 592           case (18)
 593            pop(j,18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
 594           case (19)
 595            pop(j,19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
 596          end select
 597         end do
 598        end do
 599       end subroutine GenPopulation
 600
 601 c ======================================================================
 602 c  CrossoverS subroutine
 603 c ======================================================================
 604 c  Does a simple value crossing-over on a pair
 605 c  of two given individuals, giving two new individuals.
 606 c ----------------------------------------------------------------------
 607
 608       subroutine CrossoverS(ind1,ind2)
 609        real*8 :: ind1(19),ind2(19),temp(19)
 610        integer*4 :: loc
 611        loc = 1 + int(rand(0)/(1/18.0))
 612        temp = ind2
 613        do i=(loc),19
 614         ind2(i)=ind1(i)
 615         ind1(i)=temp(i)
 616        end do
 617       end subroutine CrossoverS
 618
 619 c ======================================================================
 620 c  MutationV subroutine
 621 c ======================================================================
 622 c  Does a value mutation for the given individual (ind) with the
 623 c  given mutation ratio (mratio).
 624 c ----------------------------------------------------------------------
 625
 626       subroutine MutationV(ind,id,mratio)
 627        include 'GA.inc'
 628 c       real*8 :: rand
 629        real :: mratio
 630        real*8 :: ind(19),ti
 631        integer*4 :: i,id
 632        character*100 :: tmptext
 633        character*150 :: logtext
 634
 635        do i = 1,19
 636         if (IGNORE_WEIGHT(i).eq.1) then
 637 c            write(tmptext,'(I2)') i
 638 c           call write2log("Skipping weight "//trim(tmptext))
 639             goto 2300
 640         endif
 641         if (rand(0)<mratio) then
 642           select case(i)
 643             case (1)
 644              ti=ind(1)
 645              ind(1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
 646             case (2)
 647              ti=ind(2)
 648              ind(2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
 649             case (3)
 650              ti=ind(3)
 651              ind(3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
 652             case (4)
 653              ti=ind(4)
 654              ind(4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
 655             case (5)
 656              ti=ind(5)
 657              ind(5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
 658             case (6)
 659              ti=ind(6)
 660              ind(6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
 661             case (7)
 662              ti=ind(7)
 663              ind(7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
 664             case (8)
 665              ti=ind(8)
 666              ind(8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
 667             case (9)
 668              ti=ind(9)
 669              ind(9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
 670             case (10)
 671              ti=ind(10)
 672              ind(10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
 673             case (11)
 674              ti=ind(11)
 675              ind(11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
 676             case (12)
 677              ti=ind(12)
 678              ind(12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
 679             case (13)
 680              ti=ind(13)
 681              ind(13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
 682             case (14)
 683              ti=ind(14)
 684              ind(14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
 685             case (15)
 686              ti=ind(15)
 687              ind(15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
 688             case (16)
 689              ti=ind(16)
 690              ind(16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
 691             case (17)
 692              ti=ind(17)
 693              ind(17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
 694             case (18)
 695              ti=ind(18)
 696              ind(18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
 697             case (19)
 698              ti=ind(19)
 699              ind(19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
 700         end select
 701          write(tmptext,'(I4,A13,I2,A7,F7.5,A5,F7.5)') id," at position "&
 702      &,i," from: ",ti," to: ",ind(i)
 703          logtext = "Mutation occured in individual "//tmptext
 704          call write2log (logtext)
 705
 706         endif
 707 2300   end do
 708       end subroutine
 709
 710 c ======================================================================
 711 c  FunkcjaCeluB function
 712 c ======================================================================
 713 c  Simple (dummy) scoring function
 714 c ----------------------------------------------------------------------
 715
 716       real*8 function FunkcjaCeluB(osobnik)
 717        include 'TEST.inc'
 718        real*8 :: osobnik(20)
 719        integer*4 :: i
 720        FunkcjaCeluB = 20 -exp(-abs(WLONG-osobnik(1)))                   &
 721      & - exp(-abs(WSCP-osobnik(2)))                                     &
 722      & - exp(-abs(WELEC-osobnik(3)))                                    &
 723      & - exp(-abs(WBOND-osobnik(4)))                                    &
 724      & - exp(-abs(WANG-osobnik(5)))                                     &
 725      & - exp(-abs(WSCLOC-osobnik(6)))                                   &
 726      & - exp(-abs(WTOR-osobnik(7)))                                     &
 727      & - exp(-abs(WTORD-osobnik(8)))                                    &
 728      & - exp(-abs(WCORRH-osobnik(9)))                                   &
 729      & - exp(-abs(WCORR4-osobnik(10)))                                  &
 730      & - exp(-abs(WCORR5-osobnik(11)))                                  &
 731      & - exp(-abs(WCORR6-osobnik(12)))                                  &
 732      & - exp(-abs(WEL_LOC-osobnik(13)))                                 &
 733      & - exp(-abs(WTURN3-osobnik(14)))                                  &
 734      & - exp(-abs(WTURN4-osobnik(15)))                                  &
 735      & - exp(-abs(WTURN6-osobnik(16)))                                  &
 736      & - exp(-abs(WVDWPP-osobnik(17)))                                  &
 737      & - exp(-abs(WHPB-osobnik(18)))                                    &
 738      & - exp(-abs(WSCCOR-osobnik(19)))
 739        return
 740       end function
 741
 742 c ======================================================================
 743 c  FunkcjaCeluB2 function
 744 c ======================================================================
 745 c  Simple (dummy) scoring function
 746 c ----------------------------------------------------------------------
 747
 748       real*8 function FunkcjaCeluB2(osobnik)
 749        include 'TEST.inc'
 750        real*8 :: osobnik(20)
 751        integer*4 :: i
 752        FunkcjaCeluB2 = 1 + 4*(WLONG-osobnik(1))**2                      &
 753      & + 4*(WSCP-osobnik(2))**2                                         &
 754      & + 4*(WELEC-osobnik(3))**2                                        &
 755      & + 4*(WBOND-osobnik(4))**2                                        &
 756      & + 4*(WANG-osobnik(5))**2                                         &
 757      & + 4*(WSCLOC-osobnik(6))**2                                       &
 758      & + 4*(WTOR-osobnik(7))**2                                         &
 759      & + 4*(WTORD-osobnik(8))**2                                        &
 760      & + 4*(WCORRH-osobnik(9))**2                                       &
 761      & + 4*(WCORR4-osobnik(10))**2                                      &
 762      & + 4*(WCORR5-osobnik(11))**2                                      &
 763      & + 4*(WCORR6-osobnik(12))**2                                      &
 764      & + 4*(WEL_LOC-osobnik(13))**2                                     &
 765      & + 4*(WTURN3-osobnik(14))**2                                      &
 766      & + 4*(WTURN4-osobnik(15))**2                                      &
 767      & + 4*(WTURN6-osobnik(16))**2                                      &
 768      & + 4*(WVDWPP-osobnik(17))**2                                      &
 769      & + 4*(WHPB-osobnik(18))**2                                        &
 770      & + 4*(WSCCOR-osobnik(19))**2
 771        return
 772       end function
 773
 774 c ======================================================================
 775 c  Weights2str subroutine
 776 c ======================================================================
 777 c  Writes weights of individual(osobnik) to string (lista)
 778 c ----------------------------------------------------------------------
 779
 780       subroutine Weights2Str(osobnik,coff,lista)
 781       real*8, dimension(19) :: osobnik
 782       real*8 :: coff
 783       character(80) :: a
 784       integer :: i
 785       character(80) :: lista
 786       dimension lista(5)
 787 c       Zerowanie tablicy
 788       do i=1,5
 789        lista(i)=''
 790       enddo
 791 c
 792       !write(*,*) coff
 793       write(a,"(A6,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)") "WLONG=",os&
 794      &obnik(1)," WSCP=",osobnik(2)," WELEC=",osobnik(3)," WBOND=",osobni&
 795      &k(4)," WANG=",osobnik(5),"            &"
 796       lista(1)=a
 797       write(a,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)") "WSCLOC=",o&
 798      &sobnik(6)," WTOR=",osobnik(7)," WTORD=",osobnik(8)," WCORRH=",osob&
 799      &nik(9),"WCORR5=",osobnik(11),"        &"
 800       lista(2)=a
 801       write(a,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)") "WCORR6=",o&
 802      &sobnik(12)," WEL_LOC=",osobnik(13)," WTURN3=",osobnik(14)," WTURN4&
 803      &=",osobnik(15),"WTURN6=",osobnik(16),"    &"
 804       lista(3)=a
 805       write(a,"(A7,F7.5,A6,F7.5,A8,F7.5,A)") "WVDWPP=",osobnik(17)," WHP&
 806      &B=",osobnik(18)," WSCCOR=",osobnik(19),"                          &
 807      &           &"
 808       lista(4)=a
 809       write(a,"(A7,F7.5,A8,F7.5)") "CUTOFF=",coff," WCORR4=",osobnik(10)
 810       lista(5)=a
 811
 812       end subroutine
 813
 814 c ======================================================================
 815 c  FunkcjaCeluR function
 816 c ======================================================================
 817 c  Real scoring function
 818 c ----------------------------------------------------------------------
 819
 820       real*8 function FunkcjaCeluR(zscore)
 821        real*8 :: zscore
 822        if (zscore.eq.0) then
 823         FunkcjaCeluR = 0
 824        else
 825         FunkcjaCeluR = 1/zscore
 826        endif
 827        return
 828       end function
 829
 830 c ======================================================================
 831 c  ReadInput subroutine
 832 c ======================================================================
 833
 834       subroutine ReadInput(status)
 835       include "io.inc"
 836       include 'common.inc'
 837       include 'GA.inc'
 838       logical :: ex
 839       integer :: stat
 840       integer*4 :: status
 841       character*100 :: wiersz,tmp
 842       inquire(FILE=inpfn,EXIST=ex)
 843       if (ex) then
 844        status = 0
 845        call write2log('')
 846        call write2log('==[ Reading main input ]======================')
 847        call write2log('')
 848        open(inp,file=inpfn)
 849        do
 850          read(inp, '(A)', iostat=stat) wiersz
 851          if (stat /= 0) exit
 852 c PDB=
 853          if ((wiersz(1:4).eq.'PDB=').or.(wiersz(1:4).eq.'pdb=')) then
 854           npdb=1
 855           tmp = wiersz(5:len_trim(wiersz))
 856           do i=1,len_trim(tmp)
 857            if (tmp(i:i).eq.' ') then
 858             npdb=npdb+1
 859            endif
 860           end do
 861           if (npdb.gt.maxnpdb) then
 862            call write2log("Number of input PDB exceeds maxnpdb!")
 863            status = 1
 864            exit
 865           endif
 866           do i=1,npdb
 867            if (index(trim(tmp),' ').gt.0) then
 868             pdbfiles(i)=tmp(1:index(trim(tmp),' '))
 869            else
 870             pdbfiles(i)=tmp(1:len_trim(tmp))
 871            endif
 872            tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
 873           end do
 874          endif ! Koniec czytania "PDB="
 875
 876 c ALG=
 877          if ((wiersz(1:4).eq.'ALG=').or.(wiersz(1:4).eq.'alg=')) then
 878           alg=wiersz(5:len_trim(wiersz))
 879           select case(alg)
 880            case('simple')
 881             call write2log ("Using 'simple' algorithm")
 882            case('better')
 883             call write2log ("Using 'better' algorithm")
 884            case('betterv2')
 885             call write2log ("Using 'betterv2' algorithm")
 886            case('csa')
 887             call write2log ("Using 'CSA' algorithm")
 888            case('cluster')
 889             call write2log ("Using 'Clustering CSA' algorithm")
 890            case default
 891             alg="csa"
 892             call write2log ("Unknown algorithm. Using 'csa' as default")
 893           end select
 894          endif
 895
 896 c GENERATIONS=
 897          select case (wiersz(1:12))
 898           case('GENERATIONS=','generations=')
 899            tmp = wiersz(13:len_trim(wiersz))
 900            read(tmp,'(I4)') maxgen
 901          end select
 902
 903 c CICUTOFF=
 904 c         select case(wiersz(1:9))
 905 c          case('CICUTOFF=','cicutoff=')
 906 c          tmp = wiersz(10:len_trim(wiersz))
 907 c          read(tmp(1:len_trim(tmp)),'(F7.5)') cicutoff
 908 c          call write2log("Initial CSA cutoff is set to "//tmp)
 909 c         end select
 910
 911 c MINCO=
 912          select case(wiersz(1:6))
 913           case('MINCO=','minco=')
 914           tmp = wiersz(7:len_trim(wiersz))
 915           read(tmp(1:len_trim(tmp)),'(F7.5)') minco
 916           call write2log("Minimal CSA cutoff factor is set to "//tmp)
 917          end select
 918
 919 c MAXCO=
 920          select case(wiersz(1:6))
 921           case('MAXCO=','maxco=')
 922           tmp = wiersz(7:len_trim(wiersz))
 923           read(tmp(1:len_trim(tmp)),'(F7.5)') maxco
 924           call write2log("Maximal CSA cutoff factor is set to "//tmp)
 925          end select
 926
 927
 928 c POPULATION=
 929          select case(wiersz(1:11))
 930           case('POPULATION=','population=')
 931            tmp = wiersz(12:len_trim(wiersz))
 932            read(tmp,'(I4)') banksize
 933            call write2log("Bank size is set to "//tmp)
 934          end select
 935
 936 c WHAMTEMPLATE
 937          select case(wiersz(1:13))
 938           case('WHAMTEMPLATE=','whamtemplate=')
 939            tmp = wiersz(14:len_trim(wiersz))
 940            ntwham=1
 941            do i=1,len_trim(tmp)
 942             if (tmp(i:i).eq.' ') then
 943              ntwham=ntwham+1
 944             endif
 945            end do
 946            if (ntwham.gt.maxnpdb) then
 947             call write2log("Number of WHAM templates exceeds maxnpdb!")
 948              status = 1
 949              exit
 950            endif
 951            do i=1,ntwham
 952             if (index(trim(tmp),' ').gt.0) then
 953               whamtemplate(i)=tmp(1:index(trim(tmp),' '))
 954             else
 955               whamtemplate(i)=tmp(1:len_trim(tmp))
 956             endif
 957               tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
 958             end do
 959          end select
 960
 961 c MREMDTEMPLATE=
 962          select case(wiersz(1:14))
 963           case('MREMDTEMPLATE=','mremdtemplate=')
 964           tmp = wiersz(15:len_trim(wiersz))
 965           ntmremd=1
 966           do i=1,len_trim(tmp)
 967            if (tmp(i:i).eq.' ') then
 968            ntmremd=ntmremd+1
 969            endif
 970           end do
 971           if (ntmremd.gt.maxnpdb) then
 972             call write2log("Number of MREMD templates exceeds maxnpdb!")
 973             status = 1
 974             exit
 975           endif
 976           do i=1,ntmremd
 977             if (index(trim(tmp),' ').gt.0) then
 978              mremdtemplate(i)=tmp(1:index(trim(tmp),' '))
 979             else
 980              mremdtemplate(i)=tmp(1:len_trim(tmp))
 981             endif
 982             tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
 983           end do
 984          end select
 985
 986 c MAXMIN=
 987          select case(wiersz(1:7))
 988           case('MAXMIN=','maxmin=')
 989           tmp = wiersz(8:len_trim(wiersz))
 990           read(tmp,'(I4)') maxminstep
 991           if (maxminstep.gt.0) then
 992            call write2log("ZSCORE weights minimalization set to "//tmp)
 993            do_optima=.true.
 994           else
 995            call write2log("ZSCORE weights minimalization disabled")
 996            do_optima=.false.
 997           endif
 998          end select
 999
1000 c SCRIPTS=
1001          select case(wiersz(1:8))
1002           case('SCRIPTS=','scripts=')
1003           nscripts=1
1004           tmp = wiersz(9:len_trim(wiersz))
1005           do i=1,len_trim(tmp)
1006            if (tmp(i:i).eq.' ') then
1007             nscripts=nscripts+1
1008            endif
1009           end do
1010           if (nscripts.gt.maxscripts) then
1011            call write2log("Number of scripts exceeds maxscripts!")
1012            status = 1
1013            exit
1014           endif
1015           do i=1,nscripts
1016            if (index(trim(tmp),' ').gt.0) then
1017             scripts(i)=tmp(1:index(trim(tmp),' '))
1018            else
1019             scripts(i)=tmp(1:len_trim(tmp))
1020            endif
1021            call write2log("Script will be used: "//scripts(i))
1022            tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
1023           end do
1024          end select ! Koniec czytania "scripts="
1025         end do
1026
1027         close(inp)
1028
1029 c write additional info to log
1030         write(tmp,'(F8.5)') MUTRATIO
1031         call write2log("Mutation ratio is set to "//tmp)
1032         write(tmp,'(I2)') BANK_MULTIPLIER
1033         call write2log("Bank multiplier is set to "//tmp)
1034         write(tmp,'(I5)') BANK_MULTIPLIER*banksize
1035         tmp=adjustl(tmp)
1036         call write2log("Trail population size will be "//trim(tmp)//" in&
1037      &dividuals")
1038       else
1039         call write2log("Input file unresga.inp does not exist!!")
1040         write(*,*) "Input file unresga.inp does not exist!!"
1041         status = 1
1042       endif
1043
1044 c Checking set parameters after reading input
1045
1046 c PDB= ?
1047       if (npdb.eq.0) then
1048         call write2log("Can't find 'PDB=' entry in input file")
1049         status = 1
1050       else
1051        do i=1,npdb
1052         inquire(FILE=trim(pdbfiles(i)),EXIST=ex)
1053         if (.not.ex) then
1054          call write2log("Can't find file "//trim(pdbfiles(i)))
1055          status = 1
1056         endif
1057        enddo
1058       endif
1059
1060 c GENERATIONS= ?
1061       if (maxgen.eq.0) then
1062        call write2log("Couldn't find GENERATIONS= entry in input file")
1063        status = 1
1064       endif
1065
1066 c POPULATION= ?
1067       if (banksize.eq.0) then
1068        call write2log("Can't find POPULATION= entry in input file")
1069        status=1
1070       elseif(mod(banksize,2).eq.1) then
1071        call write2log("POPULATION has to be a even number.")
1072        status=1
1073       endif
1074
1075 c WHAMTEMPLATE= ?
1076       if (ntwham.ne.npdb) then
1077        call write2log("Number of WHAM templates and PDB files is not equ&
1078      &al!")
1079        status=1
1080       endif
1081
1082 c MREMDTEMPLATE= ?
1083       if (ntmremd.ne.npdb) then
1084        call write2log("Number of MREMD templates and PDB files is not eq&
1085      &ual")
1086        status=1
1087       endif
1088
1089 c SCRIPTS= ?
1090       if (nscripts.gt.0) then
1091        do i=1,nscripts
1092         inquire(FILE=trim(scripts(i)),EXIST=ex)
1093         if (.not.ex) then
1094          call write2log("Can't find file "//trim(scripts(i)))
1095          status = 1
1096         endif
1097        enddo
1098       endif
1099
1100
1101       end subroutine ReadInput
1102
1103
1104 c ======================================================================
1105 c  CreateInputs subroutine
1106 c ======================================================================
1107 c  Creates input files for MREMD, WHAM and ZSCORE.
1108 c ----------------------------------------------------------------------
1109
1110       subroutine CreateInputs(rozmiar,pop)
1111        include 'io.inc'
1112        include 'GA.inc'
1113        include 'common.inc'
1114        character(80),dimension(5) :: wagi
1115        integer*4 :: rozmiar,i,j,k
1116        real*8 :: pop
1117        dimension pop(rozmiar,21)
1118        !character*50 :: pdbfiles(maxnpdb)
1119        character*50 :: command,tmptext, prefix
1120        character*256 :: line,plik
1121        integer :: gdzie = 0
1122        character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1123      &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&
1124      &.5,F20.5)"
1125        !write(tmptext,'(I10)') banksize
1126        !call write2log(trim(tmptext))
1127
1128         do j=1,npdb
1129          prefix=pdbfiles(j)(1:len_trim(pdbfiles(j))-4)
1130 c Create MREMD input
1131          open(imremd, file = mremdtemplate(j))
1132          plik=trim(prefix)//"_"//omremdfn
1133          call write2log(trim("Creating MREMD input: "//plik))
1134          open(omremd, file = plik)
1135 3000     read(imremd, '(A)', end = 3010) line
1136          gdzie = index (line, '{PREFIX}')
1137          if (gdzie .ne. 0) then
1138           line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1139           write(omremd,"(A)") trim(line)
1140          elseif (index(line,'{POPSIZE}') .ne. 0) then
1141           gdzie = index (line, '{POPSIZE}')
1142           write(tmptext,'(I4)') rozmiar
1143           line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(
1144      &gdzie+9:))
1145           write(omremd,"(A)") trim(line)
1146          elseif (index(line,'{WEIGHTS}') .ne. 0) then
1147           gdzie = index (line, '{WEIGHTS}')
1148           do i=1,rozmiar
1149            call Weights2Str(pop(i,:),CUTOFF,wagi)
1150            write(omremd,"(A)") wagi(1)
1151            write(omremd,"(A)") wagi(2)
1152            write(omremd,"(A)") wagi(3)
1153            write(omremd,"(A)") wagi(4)
1154            write(omremd,"(A)") wagi(5)
1155           enddo
1156          else
1157           write(omremd,"(A)") trim(line)
1158          endif
1159
1160          goto 3000
1161 3010     continue
1162          close(imremd)
1163          close(omremd)
1164
1165 c Create WHAM input
1166          open(iwham, file = whamtemplate(j))
1167          plik=trim(prefix)//"_"//owhamfn
1168          call write2log(trim("Creating WHAM input: "//plik))
1169          open(owham, file = plik)
1170 3015     read(iwham, '(A)', end = 3020) line
1171          gdzie = index (line, '{PREFIX}')
1172          if (gdzie .ne. 0) then
1173           line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1174           write(owham,"(A)") trim(line)
1175          elseif (index(line,'{POPSIZE}') .ne. 0) then
1176           gdzie = index (line, '{POPSIZE}')
1177           write(tmptext,'(I4)') rozmiar
1178           line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1179      &gdzie+9:))
1180           write(owham,"(A)") trim(line)
1181          elseif (index(line,'{WEIGHTS}') .ne. 0) then
1182           gdzie = index (line, '{WEIGHTS}')
1183           do i=1,rozmiar
1184            call Weights2Str(pop(i,:),CUTOFF,wagi)
1185            write(owham,"(A)") wagi(1)
1186            write(owham,"(A)") wagi(2)
1187            write(owham,"(A)") wagi(3)
1188            write(owham,"(A)") wagi(4)
1189            write(owham,"(A)") wagi(5)
1190            write(owham,"(A)") ''
1191           enddo
1192          else
1193            write(owham,"(A)") trim(line)
1194          endif
1195
1196          goto 3015
1197 3020     continue
1198          close(iwham)
1199          close(owham)
1200          enddo
1201 c Create ZSCORE input
1202          open(izscore, file = izscorefn)
1203          call write2log(trim("Creating ZSCORE input: "//ozscorefn))
1204          open(ozscore, file = ozscorefn)
1205 3025     read(izscore, '(A)', end = 3030) line
1206          gdzie = index (line, '{PREFIX')
1207          if (gdzie .ne. 0) then
1208           read(line(gdzie+7:gdzie+9),'(I2.2)') k
1209           prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1210           line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+10:)&
1211      &)
1212
1213           write(ozscore,"(A)") trim(line)
1214          elseif (index(line,'{PARAM') .ne. 0) then
1215           gdzie=index(line,'{PARAM')
1216           read(line(gdzie+6:gdzie+8),'(I2.2)') k
1217           prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1218           do i=1,rozmiar
1219            write(tmptext,'(I3.3)') i
1220            line=trim(line(1:gdzie-1))//trim(prefix)//"_wham_par"//tmptex&
1221      &t
1222            write(ozscore,"(A)") trim(line)
1223           enddo
1224          elseif (index(line,'{POPSIZE}') .ne. 0) then
1225           gdzie = index (line, '{POPSIZE}')
1226           write(tmptext,'(I4)') rozmiar
1227           line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1228      &gdzie+9:))
1229           write(ozscore,"(A)") trim(line)
1230          elseif (index(line,'{WEIGHTS}') .ne. 0) then
1231           gdzie = index (line, '{WEIGHTS}')
1232           do i=1,rozmiar
1233            call Weights2Str(pop(i,:),CUTOFF,wagi)
1234            write(ozscore,"(A)") wagi(1)
1235            write(ozscore,"(A)") wagi(2)
1236            write(ozscore,"(A)") wagi(3)
1237            write(ozscore,"(A)") wagi(4)
1238            write(ozscore,"(A)") trim(wagi(5))
1239            write(ozscore,"(A)") ''
1240           enddo
1241          elseif (index(line,'{MAXMIN}') .ne. 0) then
1242           gdzie = index (line, '{MAXMIN}')
1243           write(tmptext,'(I4)') maxminstep
1244           line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1245      &gdzie+9:))
1246           write(ozscore,"(A)") trim(line)
1247          else
1248           write(ozscore,"(A)") trim(line)
1249          endif
1250
1251          goto 3025
1252 3030     continue
1253          close(izscore)
1254          close(ozscore)
1255
1256 c Copy scripts
1257          call write2log('Copying scripts is disabled.')
1258 c        do i=1,nscripts
1259 c         ! Open output file
1260 c         plik=trim(prefix)//"/"//trim(scripts(i))
1261 c         open(oscript, file = plik)
1262 c         ! Open input file
1263 c         plik=trim(scripts(i))
1264 c         open(iscript, file = plik)
1265 c          ! rewrite from input to output
1266 c3030     read(iscript, '(A)', end = 3035) line
1267 c         write(oscript,'(A)') trim(line)
1268 c         goto 3030
1269 c3035     continue
1270 c         close(oscript)
1271 c         close(iscript)
1272 c        enddo
1273 c        enddo
1274
1275 c ==========================================
1276 c  Dump weights
1277 c ==========================================
1278
1279
1280 c       write(ow,'(I4)') generation
1281 c       write(ow,'(200A)') "# WLONG   WSCP    WELEC   WBOND    WANG   WSC&
1282 c     &LOC  WTOR    WTORD   WCORRH  WCORR4  WCORR5  WCORR6  WEL_LOC WTURN&
1283 c     &3  WTURN4  WTURN6  WVDWPP  WHPB    WSCCOR"
1284 c       do I=1,rozmiar
1285 c         write(ow,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),pop&
1286 c     &(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),pop&
1287 c     &(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,19)
1288 c        enddo
1289 c       close(ow)
1290
1291       end subroutine CreateInputs
1292
1293
1294 c ======================================================================
1295 c  ZnajdzPodobnego function
1296 c ======================================================================
1297 c  Finds the most similar individual to osobnik in population(populacja)
1298 c  and returns his position in the population
1299 c ----------------------------------------------------------------------
1300
1301       integer*4 function ZnajdzPodobnego(rozmiar,populacja,osobnik,csacu&
1302      &toff)
1303       integer*4 :: rozmiar,pozycja
1304       real*8 :: delta,mindelta,csacutoff
1305       character*30 :: tmp
1306       real*8, dimension(rozmiar,21) :: populacja
1307       real*8, dimension(21) :: osobnik
1308
1309       mindelta = csacutoff
1310 c mindelta = 10000.0
1311       pozycja = 0
1312
1313       do i=1,rozmiar
1314         delta = 0.0
1315         do j=1,19
1316           delta=delta+(populacja(i,j)-osobnik(j))**2
1317         enddo
1318 c        write(tmp,'(2F7.5 )') delta,csacutoff
1319 c        call write2log("Delta is "//tmp)
1320
1321         if (delta.lt.mindelta) then
1322           mindelta=delta
1323           pozycja=i
1324           write(tmp,'(F7.5,X,I4)') delta,pozycja
1325           call write2log("  Delta is "//tmp)
1326         endif
1327       end do
1328
1329       ZnajdzPodobnego = pozycja
1330       end function
1331
1332 c ======================================================================
1333 c  WybierzOsobnika function
1334 c ======================================================================
1335 c  Selects individual from population using roulette method
1336 c ----------------------------------------------------------------------
1337
1338       integer*4 function WybierzOsobnika(rozmiar,populacja)
1339       integer*4 :: osobnik, rozmiar
1340       real*8 :: los, partsum
1341
1342       real*8, dimension(rozmiar,21) :: populacja
1343
1344       partsum = 0.0
1345       osobnik = 0
1346       los = rand(0)
1347        do
1348          osobnik = osobnik + 1
1349          partsum = partsum + populacja(osobnik,21)
1350          if ((osobnik.eq.rozmiar).or.(partsum.ge.los)) exit
1351        end do
1352
1353        WybierzOsobnika = osobnik
1354        return
1355       end function
1356
1357 c ======================================================================
1358 c  Najlepszy function
1359 c ======================================================================
1360 c  Selects the individual with the best score from population
1361 c ----------------------------------------------------------------------
1362
1363       integer*4 function Najlepszy(rozmiar, maks, populacja)
1364       integer*4 :: rozmiar
1365       real*8 :: maks
1366       real*8, dimension(rozmiar,21) :: populacja
1367        Najlepszy = 0
1368       do i=1,rozmiar
1369        if (populacja(i,20).eq.maks) then
1370         Najlepszy = i
1371        endif
1372       end do
1373        return
1374       end function
1375
1376 c ======================================================================
1377 c  WriteState subroutine
1378 c ======================================================================
1379 c  Writes current state of calculations in case jobs have to be
1380 c  restartet
1381 c ----------------------------------------------------------------------
1382
1383       subroutine WriteState()
1384        include 'GA.inc'
1385        include 'io.inc'
1386        include 'common.inc'
1387
1388        open(ostate,file=ostatefn)
1389        write(ostate,'(I4)') generation
1390        write(ostate,'(L2)') do_optima
1391        write(ostate,'(L2)') do_ga
1392        write(ostate,'(F7.5)') avrd
1393        close(ostate)
1394       end subroutine WriteState
1395
1396 c ======================================================================
1397 c  ReadPopulation subroutine
1398 c ======================================================================
1399 c  Reads population from population.summary file
1400 c ----------------------------------------------------------------------
1401
1402       subroutine ReadPopulation(nind,pop)
1403        integer*4 :: i,nind
1404        real*8 :: pop
1405        character*8 :: dummy
1406        dimension pop(nind,21)
1407        include 'GA.inc'
1408        include 'io.inc'
1409        character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1410      &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&
1411      &.5,F20.5)"
1412        read(opopsum,'(A)') dummy
1413        do i=1,nind
1414        read(opopsum,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),p&
1415      &op(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),p&
1416      &op(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,1&
1417      &9),pop(i,20)
1418        end do
1419       end subroutine ReadPopulation
1420
1421 c ======================================================================
1422 c  ReadOptimaW subroutine
1423 c ======================================================================
1424 c  Reads weights from zscore generated weights_opt.zscorepar_XXX files
1425 c
1426 c  nind - number of individuals/files to read
1427 c  pop - array where weights are stored
1428 c ----------------------------------------------------------------------
1429
1430       subroutine ReadOptimaW(nind,pop)
1431        integer*4 :: i,nind,stat
1432
1433        real*8 :: pop
1434        character*40 :: filename
1435        character*200 :: tmptxt
1436        dimension pop(nind,21)
1437        logical :: err
1438        include 'GA.inc'
1439        include 'io.inc'
1440
1441 c Cosmetic stuff for pretty log
1442        call write2log("==[ Reading weights from zscore files ]=======")
1443        call write2log("# WLONG    WSCP     WELEC    WBOND    WANG     WS&
1444      &CLOC   WTOR     WTORD    WCORRH   WCORR4   WCORR5   WCORR6   WEL_L&
1445      &OC  WTURN3   WTURN4   WTURN6   WVDWPP   WHPB     WSCCOR")
1446
1447 c Actual weight read and error handling
1448        do i=1,nind
1449         err=.false.
1450         write(tmptxt,'(I3.3)') i
1451         filename = "weights_opt.zscorepar_"//trim(tmptxt)
1452         open(izsoptw, file = filename)
1453         read(izsoptw,"(A4,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)",iosta&
1454      &t=stat) tmptxt,pop(i,1),tmptxt,pop(i,2),tmptxt,pop(i,3),tmptxt,pop&
1455      &(i,4),tmptxt,pop(i,5)
1456         if (stat/=0) then
1457          err=.true.
1458          call write2log("-- ERROR while reading "//trim(filename)//"--")
1459          do j=1,21
1460           pop(i,j)=0.0
1461          enddo
1462          goto 4000
1463         endif
1464         read(izsoptw,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1465      &t=stat) tmptxt,pop(i,6),tmptxt,pop(i,7),tmptxt,pop(i,8),tmptxt,pop&
1466      &(i,9),tmptxt,pop(i,11)
1467         if (stat/=0) then
1468          err=.true.
1469          call write2log("-- ERROR while reading "//trim(filename)//"--")
1470          do j=1,21
1471           pop(i,j)=0.0
1472          enddo
1473          goto 4000
1474         endif
1475         read(izsoptw,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1476      &t=stat) tmptxt,pop(i,12),tmptxt,pop(i,13),tmptxt,pop(i,14),tmptxt,&
1477      &pop(i,15),tmptxt,pop(i,16)
1478          if (stat/=0) then
1479          err=.true.
1480          call write2log("-- ERROR while reading "//trim(filename)//"--")
1481          do j=1,21
1482           pop(i,j)=0.0
1483          enddo
1484          goto 4000
1485         endif
1486         read(izsoptw,"(A7,F7.5,A8,F7.5,A6,F7.5,A)",iostat=stat) tmptxt, &
1487      &pop(i,19),tmptxt,pop(i,17),tmptxt,pop(i,18),tmptxt
1488          if (stat/=0) then
1489          err=.true.
1490          call write2log("-- ERROR while reading "//trim(filename)//"--")
1491          do j=1,21
1492           pop(i,j)=0.0
1493          enddo
1494          goto 4000
1495         endif
1496         read(izsoptw,"(A22,F7.5)",iostat=stat) tmptxt, pop(i,10)
1497         if (stat/=0) then
1498         err=.true.
1499         call write2log("-- ERROR while reading "//trim(filename)//"--")
1500         do j=1,21
1501          pop(i,j)=0.0
1502         enddo
1503         endif
1504
1505         if (.not.err) then
1506          write(tmptxt,'(19F9.5)') (pop(i,j),j=1,19)
1507          call write2log(trim(tmptxt))
1508         endif
1509
1510 4000    close(izsoptw)
1511        end do
1512       end subroutine ReadOptimaW
1513
1514
1515 c ======================================================================
1516 c  ReadZEnergy subroutine
1517 c ======================================================================
1518 c  Reads Energy from ZSCORE output files zscore_out_parXXX_GB000
1519 c ----------------------------------------------------------------------
1520
1521       subroutine ReadZEnergy(nind,pop)
1522        integer*4 :: i,nind
1523        integer*4 :: stat
1524        real*8 :: pop
1525        character*40 :: filename,tmptxt,tmptxt2
1526        character*100 :: wiersz
1527        dimension pop(nind,21)
1528        include 'io.inc'
1529
1530        call write2log("==[ Reading Final Function Value from zscore ]===&
1531      &====")
1532        do i=1,nind
1533         write(tmptxt,'(I3.3)') i
1534         pop(i,20)=0
1535         filename=trim(izenergyfn)//trim(tmptxt)//"_GB000"
1536         open(izenergy, file = filename)
1537         do
1538          read(izenergy, '(A)', iostat=stat) wiersz
1539          if (stat /= 0) exit
1540          if (wiersz(1:22).eq.' Final function value:') then
1541           tmptxt=trim(adjustl(wiersz(23:48)))
1542           read(tmptxt,'(F16.5)') pop(i,20)
1543           write(tmptxt2,'(E15.7)') pop(i,20)
1544           call write2log("Reading score from "//filename//" : "//tmptxt/&
1545      &/" : "//tmptxt2)
1546          endif
1547         end do
1548         if (pop(i,20).eq.0) then
1549 c Setting FFV to a realy big valiue or otherwise CalcFittness will have problems
1550          pop(i,20)=huge(0.0d0)
1551          call write2log("ERROR while reading FFV from "//filename)
1552         endif
1553         close(izenergy)
1554        end do
1555
1556       end subroutine ReadZEnergy
1557
1558 c =======================================================================
1559 c  WritePopSum subroutine
1560 c =======================================================================
1561 c  Writes whole individuals (weights, score, ID, generation) to file
1562 c  population.summary  This file can be used to plot evolution of the
1563 c  weights and scores over time.
1564 c -----------------------------------------------------------------------
1565       subroutine WritePopSum()
1566        logical :: jestplik = .false.
1567        character*5 :: tmptext
1568 c       integer*4 :: nind
1569 c       real*8 :: bank
1570 c       dimension bank(nind,21)
1571        include 'io.inc'
1572        include 'common.inc'
1573
1574
1575 c      write(*,*)  banksize
1576
1577       write(tmptext,'(I3.3)') generation
1578       inquire(file=opopsumfn, exist=jestplik)
1579       if (.not. jestplik) then
1580        open(opopsum, file = opopsumfn)
1581        write(opopsum, "(A189)") "# WLONG   WSCP    WELEC   WBOND   WANG &
1582      &   WSCLOC  WTOR    WTORD   WCORRH  WCORR4  WCORR5  WCORR6  WEL_LOC&
1583      & WTURN3  WTURN4  WTURN6  WVDWPP  WHPB    WSCCOR                SCO&
1584      &RE    IND      GEN"
1585        close(opopsum)
1586       endif
1587
1588       open(opopsum, file = opopsumfn, access = "append")
1589 c       call ReadPopulation(banksize,populacja)
1590       write(opopsum, "(A189)") "# ==== Generation "//trim(tmptext)//" ==&
1591      &==================================================================&
1592      &==================================================================&
1593      &================================="
1594 c       call write2log(banksize)
1595        do i=1,banksize
1596 c      write(*,*) i
1597 c        write(opopsum,"(F9.5,2F8.5)") populacja(i,1),populacja(i,2),popu&
1598 c     &lacja(i,3)
1599       write(opopsum,"(F9.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,&
1600      &F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F20.5,I8,I8)") bank(i&
1601      &,1),bank(i,2),bank(i,3),bank(i,4),bank(i,5),bank(i,6),bank(i,7),ba&
1602      &nk(i,8),bank(i,9),bank(i,10),bank(i,11),bank(i,12),bank(i,13),bank&
1603      &(i,14),bank(i,15),bank(i,16),bank(i,17),bank(i,18),bank(i,19),bank&
1604      &(i,20),i,iter
1605        enddo
1606       write(opopsum,*) ""
1607       close(opopsum)
1608
1609       end subroutine WritePopSum
1610
1611 c ======================================================================
1612 c  WriteBank subroutine
1613 c ======================================================================
1614 c  Writes CSA BANK to file
1615 c ----------------------------------------------------------------------
1616       subroutine WriteBank(b)
1617        include 'io.inc'
1618        include 'common.inc'
1619        real*8,dimension(banksize,21) :: b
1620        character*250 :: tmptext
1621        character*250 :: header
1622
1623
1624        header="#    WLONG      WSCP     WELEC     WBOND      WANG    WSC&
1625      &LOC      WTOR     WTORD    WCORRH    WCORR4    WCORR5    WCORR6   &
1626      &WEL_LOC    WTURN3    WTURN4    WTURN6    WVDWPP      WHPB    WSCCO&
1627      &R            FFV   FITNESS"
1628
1629        call write2log("Writing Bank to file "//iobankfn)
1630        open(iobank, file = iobankfn)
1631        write(iobank, "(A)") trim(adjustl(header))
1632
1633        do i=1,banksize
1634         write(iobank,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1635         write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1636         call write2log(trim(tmptext))
1637        enddo
1638        close(iobank)
1639       end subroutine
1640
1641 c ======================================================================
1642 c  ReadBank subroutine
1643 c ======================================================================
1644 c  Read CSA BANK from file
1645 c ----------------------------------------------------------------------
1646       subroutine ReadBank(b)
1647        include 'io.inc'
1648        include 'common.inc'
1649        real*8,dimension(banksize,21) :: b
1650        character*250 :: tmptext
1651
1652        call write2log("Reading Bank from file "//iobankfn)
1653        open(iobank, file = iobankfn)
1654 c skip header
1655        read(iobank,"(A)") tmptext
1656        do i=1,banksize
1657         read(iobank,"(19F10.5,E15.7,F10.5)") (b(i,j),j=1,21)
1658         write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1659         call write2log(trim(tmptext))
1660        enddo
1661        close(iobank)
1662       end subroutine
1663
1664 c ======================================================================
1665 c  CalcFitness subroutine
1666 c ======================================================================
1667 c  Calculates fitness of every individual in the bank
1668 c ----------------------------------------------------------------------
1669       subroutine CalcFitness(b,fitn)
1670        include 'common.inc'
1671        real*8,dimension(banksize,21) :: b
1672        real*8 :: fitn, sumfitn
1673
1674        fitn = 0.0
1675        sumfitn = 0.0
1676
1677        do i=1,banksize
1678         fitn=fitn+(1/b(i,20))
1679         sumfitn=sumfitn+b(i,20)
1680        end do
1681        do i=1,banksize
1682         b(i,21)=1/(b(i,20)*fitn)
1683 c        b(i,21)=b(i,20)/fitn
1684        end do
1685        fitn=sumfitn
1686       end subroutine
1687
1688 c ======================================================================
1689 c  CalcAvgDist subroutine
1690 c ======================================================================
1691 c  Calculates average distance between individuals in the bank
1692 c ----------------------------------------------------------------------
1693       subroutine CalcAvgDist(b,avgd)
1694        include 'common.inc'
1695        real*8,dimension(banksize,21) :: b
1696        real*8 :: d,avgd
1697        integer*4 :: nd,w
1698
1699        d=0.0                          ! distance
1700        nd = (banksize-1)*banksize/2   ! number of distances to calculate
1701
1702        do i=1,banksize-1
1703         do j=i+1,banksize
1704          do w=1,19
1705           d=d+(b(i,w)-b(j,w))**2
1706          end do
1707         end do
1708        end do
1709        avgd=sqrt(d)/nd
1710       end subroutine
1711
1712 c -----------------------------------------------------------------------
1713
1714