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