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