website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/randomize/randomize.f

   1 c===========================================================================
   2 c
   3 c   This file is part of TISEAN
   4 c
   5 c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   6 c
   7 c   TISEAN is free software; you can redistribute it and/or modify
   8 c   it under the terms of the GNU General Public License as published by
   9 c   the Free Software Foundation; either version 2 of the License, or
  10 c   (at your option) any later version.
  11 c
  12 c   TISEAN is distributed in the hope that it will be useful,
  13 c   but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 c   GNU General Public License for more details.
  16 c
  17 c   You should have received a copy of the GNU General Public License
  18 c   along with TISEAN; if not, write to the Free Software
  19 c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 c
  21 c===========================================================================
  22 c   constrained randomization
  23 c   author T. Schreiber (1999)
  24 c===========================================================================
  25       parameter(nx=100000,mx=20)
  26       double precision time
  27       dimension x(nx,mx), y(nx,mx), xx(nx,mx), icol(mx)
  28       character*72 file, fout, comment
  29       common nmax,cost,temp,cmin,rate,x
  30       external rand
  31       data wr/0.9/, nsur/1/
  32       data iverb/15/
  33
  34       call whatido("constrained randomization",iverb)
  35       call what_cost()
  36       call what_cool()
  37       call what_permute()
  38       rr=rand(ican("I",0)/real(2**22))
  39       nsur=min(999,ican("n",nsur))
  40       nmax=ican("l",nx)
  41       nexcl=ican("x",0)
  42       mcmax=ican("m",0)
  43       call columns(mc,mx,icol)
  44       if(mcmax.eq.0) mcmax=max(1,mc)
  45       wr=fcan("u",wr)
  46       call opts_cost(mcmax)
  47       call opts_cool()
  48       call opts_permute()
  49       isout=igetout(fout,iverb)
  50
  51       call nthstring(1,file)
  52       call xreadfile(nmax,mcmax,nx,xx,nexcl,icol,file,iverb)
  53       call cost_transform(nmax,mcmax,nx,xx)
  54       if(file.eq."-") file="stdin"
  55       if(isout.eq.1) call addsuff(fout,file,"_rnd")
  56       if(nsur.gt.1) call suffix(fout,"_000")
  57
  58       do 10 isur=1,nsur
  59          do 20 n=1,nmax
  60             do 20 m=1,mcmax
  61  20            x(n,m)=xx(n,m)
  62          rate=1
  63          cost=r1mach(2)
  64          cmin=cost
  65          if(nsur.gt.1) write(fout(index(fout," ")-3:72),'(i3.3)') isur
  66          temp=cool_init()
  67          call cost_init()
  68          call permute_init()
  69          call cost_full(iv_vcost(iverb))
  70          cmin=cost
  71          time=0
  72  1       time=time+1.
  73          call permute(n1,n2)
  74          cmax=cost-temp*log(rand(0.0))    ! maximal acceptable cost
  75          call cost_update(n1,n2,cmax,iaccept,iv_vmatch(iverb))
  76          tnew=cool(iaccept,iend,iv_cool(iverb))
  77          if(tnew.ne.temp.or.cost.lt.cmin*wr) then
  78             cc=cost
  79             call cost_full(iv_vcost(iverb))
  80             if(iv_match(iverb).eq.1) write(istderr(),*)
  81      .         "cost function mismatch: ", abs((cc-cost)/cost)
  82          endif
  83          temp=tnew
  84          if(cost.lt.cmin*wr) then
  85             if(iv_cost(iverb).eq.1) write(istderr(),*)
  86      .         "after ",real(time)," steps at T=",temp," cost: ",cost
  87             cmin=cost
  88             call cost_inverse(nmax,mcmax,nx,x,y)
  89             write(comment,'(8h# cost: ,g15.5)') cost
  90             call xwritecfile(nmax,mcmax,nx,y,fout,iverb,comment)
  91          endif
  92          if(iend.ne.1) goto 1
  93          write(comment,'(8h# cost: ,g15.5)') cost
  94          call writecfile(nmax,mcmax,nx,y,fout,iverb,comment)
  95  10      continue
  96       end
  97
  98       subroutine usage()
  99 c usage message
 100
 101       call whatineed("[-n# -u# -I# -o outfile -l# -x# -c# -V# -h]"//
 102      .   " [cost opt.] [cooling opt.] [permutation opt.] file")
 103       call popt("n","number of surrogates (1)")
 104       call popt("u","improvement factor before write (0.9)")
 105       call popt("I","seed for random numbers (0)")
 106       call popt("l","number of values to be read (all)")
 107       call popt("x","number of values to be skipped (0)")
 108       call popt("c","column to be read (1 or file,#)")
 109       call pout("file_rnd(_nnn)")
 110       call pall()
 111       call ptext("Verbosity levels (add what you want):")
 112       call ptext("          1 = input/output" )
 113       call ptext("          2 = current cost if improved")
 114       call ptext("          4 = cost mismatch")
 115       call ptext("          8 = temperature etc. at cooling")
 116       call ptext("         16 = verbose cost if improved")
 117       call ptext("         32 = verbose cost mismatch")
 118       write(istderr(),'()')
 119       call usage_cost()
 120       write(istderr(),'()')
 121       call usage_cool()
 122       write(istderr(),'()')
 123       call usage_permute()
 124       write(istderr(),'()')
 125       stop
 126       end
 127