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