source/ga/GA.f

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