source/ga/GA.f

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