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