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