← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/Fortran/schmm.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/schmm.f
 
1
*----------------------------------------------------------------------*     
 
2
* $Id: schmm.f 92 2003-03-24 14:47:54Z vflegel $
 
3
*----------------------------------------------------------------------*     
 
4
*       Version:  File under developpment for release 2.3
 
5
*----------------------------------------------------------------------*     
 
6
      Subroutine SCHMM 
 
7
     *   (NOUT,IDMP,RIHM,RMHM,LHMM,NABC,FLOW,FSCA,OPTF,OPFF,RD,RI)
 
8
      
 
9
 
 
10
      Real             FNOR(26)
 
11
 
 
12
      Include         'pfind.f'
 
13
      Include         'hmdat.f'
 
14
      Include         'sterr.f'
 
15
 
 
16
      Logical          OPTF
 
17
      Logical          OPFF
 
18
 
 
19
* rescale emission probabilities
 
20
 
 
21
      Do  20 I1=0,LHMM
 
22
         Do I2=1,NABC
 
23
            FNOR(I2)=RIHM(I2,I1)
 
24
         End do 
 
25
         Call GetNP(FNOR,NABC,RNOR,FLOW)
 
26
         If(RNOR.NE.0.0) then 
 
27
            Do I2=1,NABC
 
28
               If(RIHM(I2,I1).GT.FLOW) RIHM(I2,I1)=RIHM(I2,I1)+RNOR
 
29
            End do 
 
30
            If(RIHM(IM,I1).GT.FLOW) RIHM(IM,I1)=RIHM(IM,I1)-RNOR
 
31
            If(RIHM(II,I1).GT.FLOW) RIHM(II,I1)=RIHM(II,I1)-RNOR
 
32
            If(RIHM(ID,I1).GT.FLOW) RIHM(ID,I1)=RIHM(ID,I1)-RNOR
 
33
         End if 
 
34
 
 
35
         If(RIHM(II,I1).GE.0) then 
 
36
            Write(NOUT,'(''# Insert position '',I4,
 
37
     *         '': Log Prob(I->I)='',F8.4,'' reset to -0.01'')')
 
38
     *         I1,RIHM(II,I1)
 
39
            RIHM(II,I1)=-0.01
 
40
         End if 
 
41
 
 
42
         If(I1.EQ.0) then
 
43
            Do I2=1,NABC
 
44
               RMHM(I2,I1)=FLOW
 
45
            End do
 
46
            Go to  20 
 
47
         End if
 
48
 
 
49
         Do I2=1,NABC
 
50
            FNOR(I2)=RMHM(I2,I1)
 
51
         End do 
 
52
         Call GetNP(FNOR,NABC,RNOR,FLOW)
 
53
         If(RNOR.NE.0.0) then 
 
54
            Do I2=1,NABC
 
55
               If(RMHM(I2,I1).GT.FLOW) RMHM(I2,I1)=RMHM(I2,I1)+RNOR
 
56
            End do 
 
57
            If(RIHM(MM,I1).GT.FLOW) RIHM(MM,I1)=RIHM(MM,I1)-RNOR
 
58
            If(RIHM(MI,I1).GT.FLOW) RIHM(MI,I1)=RIHM(MI,I1)-RNOR
 
59
            If(RIHM(MD,I1).GT.FLOW) RIHM(MD,I1)=RIHM(MD,I1)-RNOR
 
60
         End if
 
61
 
 
62
         If(RMHM( D,I1).GT.FLOW) then
 
63
            RNOR=-RMHM( D,I1)
 
64
            RMHM( D,I1)=0.0
 
65
            RIHM(DM,I1)=RIHM(DM,I1)-RNOR
 
66
            RIHM(DI,I1)=RIHM(DI,I1)-RNOR
 
67
            RIHM(DD,I1)=RIHM(DD,I1)-RNOR
 
68
         End if 
 
69
 
 
70
 20   Continue
 
71
 
 
72
* transitions 
 
73
 
 
74
      Do  50 I1=LHMM,0,-1
 
75
 
 
76
* - insert positions
 
77
 
 
78
         FNOR( 1)=RIHM(IM,I1)
 
79
         FNOR( 2)=RIHM(ID,I1)
 
80
         Call GetNP(FNOR, 2,RNOR,FLOW)
 
81
         If(RIHM(II,  I1).GT.FLOW)
 
82
     *      RNOR=RNOR+LOG(1.0-EXP(RIHM(II,I1)))
 
83
         If(RIHM(IM,I1).GT.FLOW) RIHM(IM,I1)=RIHM(IM,I1)+RNOR
 
84
         If(RIHM(ID,I1).GT.FLOW) RIHM(ID,I1)=RIHM(ID,I1)+RNOR
 
85
         If(RIHM(MI,I1).GT.FLOW) RIHM(MI,I1)=RIHM(MI,I1)-RNOR
 
86
         If(RIHM(DI,I1).GT.FLOW) RIHM(DI,I1)=RIHM(DI,I1)-RNOR
 
87
 
 
88
         If(I1.EQ.0) go to 50
 
89
 
 
90
* -match positions 
 
91
 
 
92
         FNOR( 1)=RIHM(MM,I1)
 
93
         FNOR( 2)=RIHM(MI,I1)
 
94
         FNOR( 3)=RIHM(MD,I1)
 
95
         Call GetNP(FNOR, 3,RNOR,FLOW)
 
96
         If(RIHM(MM,I1).GT.FLOW) RIHM(MM,I1)=RIHM(MM,I1)+RNOR
 
97
         If(RIHM(MI,I1).GT.FLOW) RIHM(MI,I1)=RIHM(MI,I1)+RNOR
 
98
         If(RIHM(MD,I1).GT.FLOW) RIHM(MD,I1)=RIHM(MD,I1)+RNOR
 
99
         J1=I1-1
 
100
         If(RIHM(MM,J1).GT.FLOW) RIHM(MM,J1)=RIHM(MM,J1)-RNOR
 
101
         If(RIHM(IM,J1).GT.FLOW) RIHM(IM,J1)=RIHM(IM,J1)-RNOR
 
102
         If(RIHM(DM,J1).GT.FLOW) RIHM(DM,J1)=RIHM(DM,J1)-RNOR
 
103
 
 
104
* - delete positions 
 
105
 
 
106
         FNOR( 1)=RIHM(DM,I1)
 
107
         FNOR( 2)=RIHM(DI,I1)
 
108
         FNOR( 3)=RIHM(DD,I1)
 
109
         Call GetNP(FNOR, 3,RNOR,FLOW)
 
110
         If(RIHM(DM,I1).GT.FLOW) RIHM(DM,I1)=RIHM(DM,I1)+RNOR
 
111
         If(RIHM(DI,I1).GT.FLOW) RIHM(DI,I1)=RIHM(DI,I1)+RNOR
 
112
         If(RIHM(DD,I1).GT.FLOW) RIHM(DD,I1)=RIHM(DD,I1)+RNOR
 
113
         J1=I1-1
 
114
         If(RIHM(MD,J1).GT.FLOW) RIHM(MD,J1)=RIHM(MD,J1)-RNOR
 
115
         If(RIHM(ID,J1).GT.FLOW) RIHM(ID,J1)=RIHM(ID,J1)-RNOR
 
116
         If(RIHM(DD,J1).GT.FLOW) RIHM(DD,J1)=RIHM(DD,J1)-RNOR
 
117
 
 
118
 50   Continue 
 
119
 
 
120
      FNOR( 1)=RIHM(MM, 0)
 
121
      FNOR( 2)=RIHM(MI, 0)
 
122
      FNOR( 3)=RIHM(MD, 0)
 
123
      Call GetNP(FNOR, 3,RNOR,FLOW)
 
124
      If(RIHM(MM, 0).GT.FLOW) RIHM(MM, 0)=RIHM(MM, 0)+RNOR
 
125
      If(RIHM(MI, 0).GT.FLOW) RIHM(MI, 0)=RIHM(MI, 0)+RNOR
 
126
      If(RIHM(MD, 0).GT.FLOW) RIHM(MD, 0)=RIHM(MD, 0)+RNOR
 
127
 
 
128
      FSCA=-RNOR
 
129
 
 
130
* add FIMs
 
131
 
 
132
      If(OPTF.OR.OPFF) then
 
133
 
 
134
         If(RD.GT.0.99999) RD=0.99999
 
135
         If(RI.GT.0.99999) RI=0.99999
 
136
 
 
137
         If(RD.GT.0.0) then
 
138
            RD=LOG(RD)
 
139
         Else
 
140
            RD=FLOW
 
141
         End if 
 
142
         If(RI.GT.0.0) then
 
143
            RI=LOG(RI)
 
144
         Else
 
145
            RI=FLOW
 
146
         End if 
 
147
 
 
148
* - state M(   0)
 
149
 
 
150
         RIHM(MD,   0)=RD
 
151
         FNOR( 1)=RIHM(MI,   0)
 
152
         FNOR( 2)=RIHM(MM,   0)
 
153
         Call GetNP(FNOR, 2,RNOR,FLOW)
 
154
         If(RIHM(MD,   0).GT.FLOW)
 
155
     *      RNOR=RNOR+LOG(1.0-EXP(RIHM(MD,   0)))
 
156
         RIHM(MI,   0)=RIHM(MI,   0)+RNOR
 
157
         RIHM(MM,   0)=RIHM(MM,   0)+RNOR
 
158
 
 
159
         RIHM(MI,   0)=RD
 
160
         FNOR( 1)=RIHM(MD,   0)
 
161
         FNOR( 2)=RIHM(MM,   0)
 
162
         Call GetNP(FNOR, 2,RNOR,FLOW)
 
163
         If(RIHM(MI,   0).GT.FLOW)
 
164
     *      RNOR=RNOR+LOG(1.0-EXP(RIHM(MI,   0)))
 
165
         RIHM(MD,   0)=RIHM(MD,   0)+RNOR
 
166
         RIHM(MM,   0)=RIHM(MM,   0)+RNOR
 
167
 
 
168
* - state I(   0)
 
169
 
 
170
         RIHM(ID,   0)=RD
 
171
         FNOR( 1)=RIHM(II,   0)
 
172
         FNOR( 2)=RIHM(IM,   0)
 
173
         Call GetNP(FNOR, 2,RNOR,FLOW)
 
174
         If(RIHM(ID,   0).GT.FLOW)
 
175
     *      RNOR=RNOR+LOG(1.0-EXP(RIHM(ID,   0)))
 
176
         RIHM(II,   0)=RIHM(II,   0)+RNOR
 
177
         RIHM(IM,   0)=RIHM(IM,   0)+RNOR
 
178
 
 
179
         RIHM(II,   0)=RI
 
180
         FNOR( 1)=RIHM(ID,   0)
 
181
         FNOR( 2)=RIHM(IM,   0)
 
182
         Call GetNP(FNOR, 2,RNOR,FLOW)
 
183
         If(RIHM(II,   0).GT.FLOW)
 
184
     *      RNOR=RNOR+LOG(1.0-EXP(RIHM(II,   0)))
 
185
         RIHM(ID,   0)=RIHM(ID,   0)+RNOR
 
186
         RIHM(IM,   0)=RIHM(IM,   0)+RNOR
 
187
 
 
188
* - state M(LHMM)
 
189
 
 
190
         RIHM(ID,LHMM)=RD
 
191
         FNOR( 1)=RIHM(II,LHMM)
 
192
         FNOR( 2)=RIHM(IM,LHMM)
 
193
         Call GetNP(FNOR, 2,RNOR,FLOW)
 
194
         If(RIHM(ID,LHMM).GT.FLOW)
 
195
     *      RNOR=RNOR+LOG(1.0-EXP(RIHM(ID,LHMM)))
 
196
         RIHM(II,LHMM)=RIHM(II,LHMM)+RNOR
 
197
         RIHM(IM,LHMM)=RIHM(IM,LHMM)+RNOR
 
198
 
 
199
         RIHM(II,LHMM)=RI
 
200
         FNOR( 1)=RIHM(ID,LHMM)
 
201
         FNOR( 2)=RIHM(IM,LHMM)
 
202
         Call GetNP(FNOR, 2,RNOR,FLOW)
 
203
         If(RIHM(II,LHMM).GT.FLOW)
 
204
     *      RNOR=RNOR+LOG(1.0-EXP(RIHM(II,LHMM)))
 
205
         RIHM(ID,LHMM)=RIHM(ID,LHMM)+RNOR
 
206
         RIHM(IM,LHMM)=RIHM(IM,LHMM)+RNOR
 
207
 
 
208
* - state D(LHMM)
 
209
 
 
210
         RIHM(DI,LHMM)=RI
 
211
         If(RI.GT.FLOW) then
 
212
            RIHM(DM,LHMM)=1-EXP(RI) 
 
213
         Else
 
214
            RIHM(DM,LHMM)=1.0
 
215
         End if
 
216
         
 
217
 
 
218
* - state I(LHMM)
 
219
 
 
220
         RIHM(II,LHMM)=RI
 
221
         If(RI.GT.FLOW) then
 
222
            RIHM(IM,LHMM)=1-EXP(RI) 
 
223
         Else
 
224
            RIHM(IM,LHMM)=1.0
 
225
         End if
 
226
 
 
227
      End if
 
228
 
 
229
* - internal delete states
 
230
 
 
231
      If(OPFF.AND.RD.GT.FLOW) then
 
232
         RD=LOG(EXP(RD)/(1-EXP(RD)))
 
233
         Do I1=1,LHMM-1
 
234
            RIHM(DD,I1)=RD
 
235
            FNOR( 1)=RIHM(DD,I1)
 
236
            FNOR( 2)=RIHM(DI,I1)
 
237
            FNOR( 3)=RIHM(DM,I1)
 
238
            Call GetNP(FNOR, 3,RNOR,FLOW)
 
239
            RIHM(DD,I1)=RIHM(DD,I1)+RNOR
 
240
            RIHM(DI,I1)=RIHM(DI,I1)+RNOR
 
241
            RIHM(DM,I1)=RIHM(DM,I1)+RNOR
 
242
         End do
 
243
      End if 
 
244
 
 
245
 100  Return
 
246
      End
 
247
*----------------------------------------------------------------------*
 
248
      Subroutine GetNP(FNOR,INOR,RNOR,FLOW)
 
249
 
 
250
      Real              FNOR(*)
 
251
 
 
252
      RNOR=0.0
 
253
 
 
254
      RMAX=FLOW
 
255
      Do I1=1,INOR
 
256
         RMAX=MAX(RMAX,FNOR(I1))
 
257
      End do 
 
258
 
 
259
      If(RMAX.LE.FLOW) go to 100
 
260
 
 
261
      RSUM=0.0
 
262
      RCUT=MAX(FLOW,RMAX-20.0)
 
263
      Do I1=1,INOR
 
264
         If(FNOR(I1).GT.RCUT) RSUM=RSUM+EXP(FNOR(I1)-RMAX)
 
265
      End do 
 
266
      RNOR=-RMAX-LOG(RSUM) 
 
267
 
 
268
 100  Return
 
269
      End
 
270
*----------------------------------------------------------------------*