← Back to branch summary

~ubuntu-branches/debian/sid/pftools/sid

« back to all changes in this revision

Viewing changes to src/Fortran/pfscale.f

  • Committer: Package Import Robot
  • Author(s): Andreas Tille
  • Date: 2017-03-31 14:01:39 UTC
  • Revision ID: package-import@ubuntu.com-20170331140139-povxg86r196gejfa
Tags: upstream-3+dfsg
Import upstream version 3+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Fortran/pfscale.f
 
1
*       Program pfscale 
 
2
*----------------------------------------------------------------------*     
 
3
* $Id: pfscale.f 183 2003-11-28 11:53:33Z vflegel $
 
4
*----------------------------------------------------------------------*     
 
5
*       Function: Fits paramters of an extreme value distribution to a
 
6
*                 profile score distribution. Input: sorted score list.
 
7
*       Author:   Philipp Bucher
 
8
*       Contact:  pftools@isb-sib.ch
 
9
*       Version:  File under developpment for release 2.3
 
10
*----------------------------------------------------------------------*     
 
11
* DATA
 
12
*----------------------------------------------------------------------*
 
13
 
 
14
* array sizes, I/O units
 
15
 
 
16
      Include          'ardim.f' 
 
17
 
 
18
      Parameter        (NOUT=   6)
 
19
      Parameter        (NSCL=  10)
 
20
      Parameter        (NPRF=  11)    
 
21
 
 
22
 
 
23
* function return types
 
24
 
 
25
      Integer           Xblnk
 
26
      External          Xblnk
 
27
 
 
28
* Parameters and options
 
29
 
 
30
      Character*512     FSCL  
 
31
      Character*512     FPRF  
 
32
 
 
33
      Real*8            DL
 
34
      Integer           NN
 
35
      Integer           IMNB
 
36
      Real              RP
 
37
      Real              RQ
 
38
 
 
39
* Profile fields
 
40
 
 
41
      Include          'psdat.f'
 
42
      Include          'gsdat.f'
 
43
      Include          'djdat.f'
 
44
      Include          'nodat.f'
 
45
      Include          'codat.f'
 
46
      Include          'pfdat.f'
 
47
      Include          'dfdat.f'
 
48
      Include          'sterr.f'
 
49
      Include          'cvini.f'
 
50
 
 
51
* score statistics  
 
52
 
 
53
      Real              XSCO(IDMC)
 
54
      Real              XFRQ(IDMC)
 
55
      Real              XWGT(IDMC)
 
56
      Logical           LMNB
 
57
      Logical           LLLT
 
58
      Logical           LRNM
 
59
      Logical           LEOF
 
60
 
 
61
*----------------------------------------------------------------------*
 
62
* INPUT SECTION
 
63
*----------------------------------------------------------------------*
 
64
 
 
65
      IRC=0
 
66
      LRNM=.FALSE.
 
67
      LEOF=.FALSE.
 
68
 
 
69
* read command line 
 
70
 
 
71
      Call Repar(FSCL,FPRF,LLLT,DL,NN,RP,RQ,LMNB,IMNB,IRC)
 
72
      If(IRC.NE.0) then
 
73
         Write(NERR,'(/,
 
74
     *      ''pfscale 2.3 revision 4'',//
 
75
     *      ''Usage: pfscale [ -lhLMNPQ ] [ score-list | - ] '' 
 
76
     *      ''[ profile-file ] [ parameters ]'',/
 
77
     *      )')
 
78
         Write(NERR,'(
 
79
     *      ''   options:'',/,
 
80
     *      ''    -h: print usage help text.'',/
 
81
     *      ''    -l: do not impose limit on line length.'',/
 
82
     *      ''    -L<value>:'',/
 
83
     *      ''        logarithmic base of parameters (default: 10).'',/
 
84
     *      ''    -M<value>:'',/
 
85
     *      ''        profile mode number to scale.'',/
 
86
     *      ''    -N<value>:'',/
 
87
     *      ''        database size (default: 14147368).'',/
 
88
     *      ''    -P<value>:'',/
 
89
     *      ''        upper threshold of probability range (default:'',
 
90
     *      '' 0.0001).'',/
 
91
     *      ''    -Q<value>:'',/
 
92
     *      ''        lower threshold of probability range (default:'',
 
93
     *      '' 0.000001).'',/
 
94
     *      )')
 
95
         Write(NERR,'(
 
96
     *      ''   valid (but deprecated) parameters are:'',/,
 
97
     *      ''    [L=log-base]        use option -L instead'',/,
 
98
     *      ''    [M=mode-nb]         use option -M instead'',/,
 
99
     *      ''    [N=db-size]         use option -N instead'',/,
 
100
     *      ''    [P=upper-threshold] use option -P instead'',/
 
101
     *      ''    [Q=lower-threshold] use option -Q instead'',/
 
102
     *      )')
 
103
         Call Exit(IRC)
 
104
      End if
 
105
         
 
106
      If(FSCL.NE.' '.AND.FSCL.NE.'-') then
 
107
         MSCL=NSCL 
 
108
         Open(MSCL,File=FSCL,Status='OLD',Err=900)
 
109
      Else
 
110
         MSCL=5
 
111
      End if
 
112
 
 
113
      RL=1/LOG(DL)
 
114
 
 
115
      RDBS=RL*LOG(Real(NN))
 
116
      If(RQ.NE.0) then
 
117
         EMAX=-RL*(LOG(RQ))
 
118
      Else
 
119
         EMAX=100
 
120
      End if
 
121
      EMIN=-RL*(LOG(RP))
 
122
 
 
123
      Do   5 I1=1,IDMC
 
124
         Read(MSCL,*,End= 10,Err=901) XSCO(I1)
 
125
         XFRQ(I1)=RDBS-RL*LOG(I1-0.5)
 
126
         XWGT(I1)=I1*0.5
 
127
 5    Continue
 
128
 
 
129
 10   NSCO=I1-1       
 
130
      If(NSCO.LT.1) Go to 904
 
131
 
 
132
      XSM=0
 
133
      XFM=0
 
134
      XWM=0
 
135
      Do  20 I1=1,NSCO
 
136
         If(XFRQ(I1).GE.EMIN.AND.XFRQ(I1).LE.EMAX) then 
 
137
            XSM=XSM+XWGT(I1)*XSCO(I1)
 
138
            XFM=XFM+XWGT(I1)*XFRQ(I1)
 
139
            XWM=XWM+XWGT(I1)
 
140
         End if
 
141
 20   Continue
 
142
      XSM=XSM/XWM
 
143
      XFM=XFM/XWM
 
144
 
 
145
 
 
146
      XSV=0
 
147
      XFV=0
 
148
      XCO=0
 
149
      Do  30 I1=1,NSCO
 
150
         If(XFRQ(I1).GE.EMIN.AND.XFRQ(I1).LE.EMAX) then 
 
151
            XSV=XSV+XWGT(I1)*( (XSCO(I1)-XSM)**2 ) 
 
152
            XFV=XFV+XWGT(I1)*( (XFRQ(I1)-XFM)**2 )
 
153
            XCO=XCO+XWGT(I1)*(XFRQ(I1)-XFM)*(XSCO(I1)-XSM)
 
154
         End if
 
155
 30   Continue
 
156
      XSV=(XSV/XWM)**0.5
 
157
      XFV=(XFV/XWM)**0.5
 
158
      XCO=(XCO/XWM)/(XFV*XSV)
 
159
 
 
160
      XB=XCO*XFV/XSV
 
161
      XA=XFM-XB*XSM
 
162
 
 
163
      If(FPRF.NE.' ') Go to  50 
 
164
 
 
165
* Case 1: no profile input file - print list
 
166
 
 
167
      Write(NOUT,
 
168
     *   '(''# -LogP ='',F8.4,'' + '',F12.8,'' * raw-score'')') XA,XB 
 
169
      
 
170
      Write(NOUT,'(''#'')')
 
171
      Write(NOUT,'(''#   rank raw-score  -logFreq  -logProb'')')
 
172
      Write(NOUT,'(''#'')')
 
173
 
 
174
      Do  40 I1=1,NSCO
 
175
         XPRE=XA+XB*XSCO(I1) 
 
176
         Write(NOUT,'(I8,F10.2,F10.4,F10.4)')
 
177
     *      I1,XSCO(I1),XFRQ(I1),XPRE
 
178
 40   Continue
 
179
      Go to 100
 
180
 
 
181
* Case 2: Modify profile input file
 
182
 
 
183
 50   Continue 
 
184
 
 
185
      Call REPRF
 
186
     *   (NPRF,FPRF,
 
187
     *   CPID,CPAC,CPDT,CPDE,LHDR,CHDR,LFTR,CFTR,NABC,CABC,LPRF,LPCI,
 
188
     *   BLOG,FABC,P0,
 
189
     *   CDIS,JDIP,MDIS,NDIP,
 
190
     *   CNOR,JNOP,JNOR,MNOR,NNOR,NNPR,CNTX,RNOP,
 
191
     *   JCUT,MCLE,CCUT,ICUT,JCNM,RCUT,MCUT,
 
192
     *   IDMP,CHIP,IIPP,CHMP,IMPP,CHIL,IIPL,ILIP,
 
193
     *   CHID,IIPD,CHMD,IMPD,
 
194
     *   INBP,LEOF,.FALSE.,IRC)
 
195
 
 
196
      If(IRC.NE.0) go to 100
 
197
 
 
198
* Check normalization modes (Lowest priority mode will be updated)
 
199
 
 
200
      If(LMNB) then
 
201
         Do I1=1,JNOR
 
202
            If(NNOR(I1).EQ.IMNB) then
 
203
               J1=I1
 
204
               Go to 55
 
205
            End if
 
206
         End do
 
207
         Go to 903
 
208
      Else
 
209
         J1=1
 
210
         K1=NNPR(J1)
 
211
         Do I1=2,JNOR
 
212
            If(NNPR(I1).LT.K1) then
 
213
               K1=NNPR(I1)
 
214
               J1=I1
 
215
            End if 
 
216
         End do 
 
217
      End if
 
218
 
 
219
 55   CNTX(J1)='-LogE'
 
220
      If(MNOR(J1).NE.1) go to 902
 
221
 
 
222
* add normalisation parameters: 
 
223
 
 
224
      RNOP(1,J1)=XA
 
225
      RNOP(2,J1)=XB
 
226
 
 
227
* define cut-offs:
 
228
 
 
229
      Do I1=1,JCUT
 
230
         Do I2=1,JCNM(I1)
 
231
            If(MCUT(I2,I1).EQ.NNOR(J1)) then
 
232
               Call NtoR
 
233
     *            (RCUT(I2,I1),ICUT(I1),RNOP,KNPM,MAXN,J1,1,LSEQ,RAVE)
 
234
               Go to 60
 
235
            End if
 
236
         End do
 
237
 60      Continue
 
238
      End do 
 
239
 
 
240
* rescaling command
 
241
 
 
242
      LFTR=LFTR+1
 
243
      If(LFTR.GT.1024) LFTR=1024
 
244
      Do I1=LFTR,2,-1
 
245
         CFTR(I1)=CFTR(I1-1)
 
246
      End do
 
247
 
 
248
      CFTR(1)='CC   /RESCALED_BY="'
 
249
      Call Recmd(CFTR(1)(21:130))
 
250
      IC=Lblnk(CFTR(1))
 
251
      CFTR(1)(IC+1:)='";'
 
252
 
 
253
* write profile
 
254
 
 
255
      Call WRPRF
 
256
     *   (NOUT,LLLT,LRNM,
 
257
     *   CPID,CPAC,CPDT,CPDE,LHDR,CHDR,LFTR,CFTR,NABC,CABC,LPRF,LPCI,
 
258
     *   CDIS,JDIP,MDIS,NDIP,
 
259
     *   CNOR,JNOP,JNOR,MNOR,NNOR,NNPR,CNTX,RNOP,
 
260
     *   JCUT,MCLE,CCUT,ICUT,JCNM,RCUT,MCUT,
 
261
     *   IDMP,CHIP,IIPP,CHMP,IMPP,
 
262
     *   BLOG,FABC,P0,
 
263
     *   CHID,IIPD,CHMD,IMPD,
 
264
     *   IRC)
 
265
 
 
266
 100  Call Exit(IRC)
 
267
 
 
268
* errors
 
269
 
 
270
 900  Write(NERR,*) 'Error: Unable to open score file'//
 
271
     *   ' ''',FSCL(1:Lblnk(FSCL)),'''.'
 
272
      IRC=1
 
273
      Go to 100
 
274
 901  Write(NERR,*) 'Error: Unable to read score value.'
 
275
      IRC=1
 
276
      Go to 100
 
277
 902  Write(NERR,*) 'Error: Normalization mode number',NNOR(J1),
 
278
     *   ' is not linear. Scaling impossible.'
 
279
      IRC=1
 
280
      Go to 100
 
281
 903  Write(NERR,*) 'Error: Normalization mode number',IMNB,
 
282
     *   ' is not defined in profile.'
 
283
      IRC=1
 
284
      Go to 100
 
285
 904  Write(NERR,*) 'Error: Score list does not contain any ',
 
286
     *   'valid scores.'
 
287
      IRC=1
 
288
      Go to 100
 
289
 
 
290
      End
 
291
*----------------------------------------------------------------------*
 
292
      Subroutine Repar(FSCL,FPRF,LLLT,DL,NN,RP,RQ,LMNB,IMNB,IRC) 
 
293
 
 
294
      Include          'sterr.f'
 
295
 
 
296
      Character*(*)     FSCL  
 
297
      Character*(*)     FPRF
 
298
 
 
299
      Real*8            DL
 
300
      Integer           NN
 
301
      Real              RP
 
302
      Real              RQ
 
303
      Logical           LLLT
 
304
      Logical           LMNB
 
305
      Integer           IMNB
 
306
 
 
307
      Character*512     CARG 
 
308
 
 
309
* initializations
 
310
 
 
311
      LLLT=.TRUE.
 
312
      LLLT=.FALSE.
 
313
      IMNB=1
 
314
      FSCL=' '
 
315
      FPRF=' '
 
316
      DL=10.0
 
317
      NN=14147368
 
318
      RP=0.0001
 
319
      RQ=0.000001
 
320
      IRC=0
 
321
 
 
322
* interpret command line arguments 
 
323
 
 
324
      N1=Iargc()
 
325
 
 
326
      K1=0
 
327
      I2=1
 
328
      Do I1=1,N1
 
329
         Call GetArg(I2,CARG)
 
330
         If     (CARG(1:1).EQ.'-'.AND.CARG(2:2).NE.' ') then
 
331
            If(Index(CARG,'l').NE.0) LLLT=.FALSE.
 
332
            If(Index(CARG,'h').NE.0) go to 900
 
333
            If(Index(CARG,'L').NE.0) then
 
334
               If(CARG(3:3).NE.' ') then
 
335
                  If(CARG(3:3).EQ.'e'.OR.CARG(3:3).EQ.'E') then
 
336
                     DL=EXP(1.0) 
 
337
                  Else
 
338
                     Read(CARG(3:),*,Err=900) DL
 
339
                  End if
 
340
               Else
 
341
                  I2=I2+1
 
342
                  Call GetArg(I2,CARG)
 
343
                  If(CARG(1:1).EQ.'e'.OR.CARG(1:1).EQ.'E') then
 
344
                     DL=EXP(1.0) 
 
345
                  Else
 
346
                     Read(CARG,*,Err=900) DL
 
347
                  End if
 
348
               End if
 
349
            End if
 
350
            If(Index(CARG,'N').NE.0) then
 
351
               If(CARG(3:3).NE.' ') then
 
352
                  Read(CARG(3:),*,Err=900) NN
 
353
               Else
 
354
                  I2=I2+1
 
355
                  Call GetArg(I2,CARG)
 
356
                  Read(CARG,*,Err=900) NN
 
357
               End if
 
358
            End if
 
359
            If(Index(CARG,'P').NE.0) then
 
360
               If(CARG(3:3).NE.' ') then
 
361
                  Read(CARG(3:),*,Err=900) RP
 
362
               Else
 
363
                  I2=I2+1
 
364
                  Call GetArg(I2,CARG)
 
365
                  Read(CARG,*,Err=900) RP
 
366
               End if
 
367
            End if
 
368
            If(Index(CARG,'Q').NE.0) then
 
369
               If(CARG(3:3).NE.' ') then
 
370
                  Read(CARG(3:),*,Err=900) RQ
 
371
               Else
 
372
                  I2=I2+1
 
373
                  Call GetArg(I2,CARG)
 
374
                  Read(CARG,*,Err=900) RQ
 
375
               End if
 
376
            End if
 
377
            If(Index(CARG,'M').NE.0) then
 
378
               If(CARG(3:3).NE.' ') then
 
379
                  Read(CARG(3:),*,Err=900) IMNB
 
380
                  LMNB=.TRUE.
 
381
               Else
 
382
                  I2=I2+1
 
383
                  Call GetArg(I2,CARG)
 
384
                  Read(CARG,*,Err=900) IMNB
 
385
                  LMNB=.TRUE.
 
386
               End if
 
387
            End if
 
388
 
 
389
         Else if(CARG(1:3).EQ.'L=e'.OR.CARG(1:3).EQ.'L=E') then
 
390
            DL=EXP(1.0) 
 
391
         Else if(CARG(1:2).EQ.'L=') then
 
392
            Read(CARG(3:),*,Err=900) DL 
 
393
         Else if(CARG(1:2).EQ.'N=') then
 
394
            Read(CARG(3:),*,Err=900) NN 
 
395
         Else if(CARG(1:2).EQ.'P=') then
 
396
            Read(CARG(3:),*,Err=900) RP
 
397
         Else if(CARG(1:2).EQ.'Q=') then
 
398
            Read(CARG(3:),*,Err=900) RQ
 
399
         Else if(CARG(1:2).EQ.'M=') then
 
400
            Read(CARG(3:),*,Err=900) IMNB
 
401
            LMNB=.TRUE.
 
402
         Else if(K1.LE.1) then
 
403
            K1=K1+1
 
404
            If     (K1.EQ.1) then
 
405
               FSCL=CARG
 
406
            Else if(K1.EQ.2) then
 
407
               FPRF=CARG
 
408
            End if
 
409
C         Else
 
410
C            Go to 900
 
411
         End if
 
412
         I2=I2+1
 
413
         If(I2.GT.N1) Go to 20
 
414
      End do
 
415
      
 
416
 20   If(DL.LE.1) Go to 901
 
417
      If(NN.LT.1) Go to 902
 
418
      If(RP.LT.0.OR.RP.GT.1) Go to 903
 
419
      If(RQ.LT.0.OR.RQ.GT.1) Go to 904
 
420
      If(LMNB.AND.IMNB.LT.1) Go to 905
 
421
 100  Return
 
422
 900  IRC=1
 
423
      Go to 100
 
424
 901  Write(NERR,*) 'Error: Value of option -L is out of bound.'
 
425
      IRC=1
 
426
      Go to 100
 
427
 902  Write(NERR,*) 'Error: Value of option -N is out of bound.'
 
428
      IRC=1
 
429
      Go to 100
 
430
 903  Write(NERR,*) 'Error: Value of option -P is out of bound.'
 
431
      IRC=1
 
432
      Go to 100
 
433
 904  Write(NERR,*) 'Error: Value of option -Q is out of bound.'
 
434
      IRC=1
 
435
      Go to 100
 
436
 905  Write(NERR,*) 'Warning: Mode numbers should be greater or ',
 
437
     *   'equal to 1.'
 
438
      Go to 100
 
439
      End
 
440
*----------------------------------------------------------------------*
 
441
      Include          'reprf.f'
 
442
      Include          'wrprf.f'
 
443
      Include          'recmd.f'
 
444
      Include          'NtoR.f'
 
445
      Include          'lblnk.f'
 
446
      Include          'Xblnk.f'