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