← Back to branch summary

~ubuntu-branches/ubuntu/quantal/genometools/quantal-backports

« back to all changes in this revision

Viewing changes to src/ltr/myxdrop.gen

  • Committer: Package Import Robot
  • Author(s): Sascha Steinbiss
  • Date: 2012-07-09 14:10:23 UTC
  • Revision ID: package-import@ubuntu.com-20120709141023-juuu4spm6chqsf9o
Tags: upstream-1.4.1
Import upstream version 1.4.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/ltr/myxdrop.gen
 
1
/*
 
2
  Copyright (c) 2007 David Ellinghaus <d.ellinghaus@ikmb.uni-kiel.de>
 
3
  Copyright (c) 2007 Center for Bioinformatics, University of Hamburg
 
4
 
 
5
  Permission to use, copy, modify, and distribute this software for any
 
6
  purpose with or without fee is hereby granted, provided that the above
 
7
  copyright notice and this permission notice appear in all copies.
 
8
 
 
9
  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 
10
  WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 
11
  MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 
12
  ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 
13
  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 
14
  ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 
15
  OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 
16
*/
 
17
 
 
18
EVALXDROPARBITSCORES
 
19
{
 
20
  int integermax, integermin, i ,     /* index i */
 
21
    j = 0,                            /* index j */
 
22
    k = 0,                            /* diagonal */
 
23
    end_k = ulen - vlen,              /* diagonal of endpoint (ulen, vlen) */
 
24
    lbound = 0,                       /* diagonal lower bound */
 
25
    ubound = 0,                       /* diagonal upper bound */
 
26
    lboundtmp = 0,
 
27
    uboundtmp = 0,
 
28
    tmprow;                           /* temporary row index */
 
29
  GtUchar a, b;
 
30
  Myfrontvalue tmpfront;
 
31
  Arbitrarydistances arbitdistances;
 
32
  int d = 0,                         /* distance */
 
33
      d_pre = 0,                     /* distance previous */
 
34
      dback;
 
35
  Xdropscore bigt_tmp;               /* best score T' seen already */
 
36
  GtArrayXdropscore big_t;           /* array T[d] of best score T for each d */
 
37
  bool alwaysMININFINITYINT = true;
 
38
  int allowedMININFINITYINTgenerations = 0,
 
39
      currentMININFINITYINTgeneration = 0;
 
40
 
 
41
  gt_calculatedistancesfromscores(arbitscores, &arbitdistances);
 
42
  gt_calculateallowedMININFINITYINTgenerations(
 
43
     &allowedMININFINITYINTgenerations,
 
44
     &arbitdistances);
 
45
 
 
46
  dback = SETDBACK(xdropbelowscore);
 
47
 
 
48
  GT_INITARRAY (&big_t, Xdropscore);
 
49
  integermax = MAX (ulen, vlen);
 
50
  integermin = -integermax;
 
51
 
 
52
  /* "phase" 0 */
 
53
  for (i = 0; i < MIN(ulen,vlen); i++)
 
54
  {
 
55
    COMPARESYMBOLSSEP(i,i);
 
56
  }
 
57
 
 
58
  /* alignment already finished */
 
59
  if (i >= ulen || i >= vlen )
 
60
  {
 
61
    lbound =  1;
 
62
    ubound = -1;
 
63
  }
 
64
 
 
65
  tmpfront.dptabrow = i;
 
66
  tmpfront.dptabdirection = (GtUchar) 0;   /* no predecessor */
 
67
 
 
68
  STOREINARRAYFRONTS (fronts, 0, Myfrontvalue, tmpfront);
 
69
  xdropbest->score = bigt_tmp = S (i + i, 0);
 
70
  xdropbest->ivalue = xdropbest->jvalue = (unsigned int) i;
 
71
 
 
72
  GT_STOREINARRAY (&big_t, Xdropscore, 10, bigt_tmp);
 
73
 
 
74
  /* "phase" d > 0 */
 
75
  while (lbound <= ubound)
 
76
  {
 
77
    d++;
 
78
    d_pre = d - dback;
 
79
 
 
80
    /* calculate fronts */
 
81
    for (k = lbound - 1; k <= ubound + 1; k++)
 
82
    {
 
83
 
 
84
      i = MINUSINFINITYINT;
 
85
      /* case 1 : DELETION-EDGE  */
 
86
      if ( (lbound < k) &&
 
87
           (d - arbitdistances.del >= 0) &&
 
88
           (( -(d - arbitdistances.del) <= k - 1) &&
 
89
            (k - 1 <= d - arbitdistances.del)))
 
90
      {
 
91
        i =
 
92
          fronts->spaceMyfrontvalue[ACCESSTOFRONT(d - arbitdistances.del,
 
93
                                                  k - 1)].dptabrow + 1;
 
94
        tmpfront.dptabdirection = MYDELETIONBIT;
 
95
      }
 
96
      /* case 2: REPLACEMENT-EDGE */
 
97
      if ((lbound <= k) && (k <= ubound) && (d - arbitdistances.mis >= 0) &&
 
98
          (( -(d - arbitdistances.mis) <= k) && (k <= d - arbitdistances.mis)))
 
99
      {
 
100
        /* test, if case 1 has happened. */
 
101
        if (tmpfront.dptabdirection & MYDELETIONBIT)
 
102
        {
 
103
          tmprow =
 
104
            fronts->spaceMyfrontvalue[ACCESSTOFRONT(d - arbitdistances.mis,
 
105
                                                    k)].dptabrow + 1;
 
106
          if (tmprow > i)
 
107
          {
 
108
            i = tmprow;
 
109
            tmpfront.dptabdirection = MYREPLACEMENTBIT;
 
110
          }
 
111
        }
 
112
        else
 
113
        {
 
114
          i =
 
115
            fronts->spaceMyfrontvalue[ACCESSTOFRONT(d - arbitdistances.mis, k)]
 
116
            .dptabrow + 1;
 
117
          tmpfront.dptabdirection = MYREPLACEMENTBIT;
 
118
        }
 
119
 
 
120
      }
 
121
      /* case 3: INSERTION-EDGE */
 
122
      if ((k < ubound) && (d - arbitdistances.ins >= 0) &&
 
123
          (( -(d - arbitdistances.ins) <= k + 1) &&
 
124
           (k + 1 <= d - arbitdistances.ins)))
 
125
      {
 
126
        if ((tmpfront.dptabdirection & MYDELETIONBIT)
 
127
            || (tmpfront.dptabdirection & MYREPLACEMENTBIT))
 
128
        {
 
129
          tmprow =
 
130
            fronts->spaceMyfrontvalue[ACCESSTOFRONT (d - arbitdistances.ins,
 
131
                                                     k + 1)].dptabrow;
 
132
          if (tmprow > i)
 
133
          {
 
134
            i = tmprow;
 
135
            tmpfront.dptabdirection = MYINSERTIONBIT;
 
136
          }
 
137
        }
 
138
        else
 
139
        {
 
140
          i =
 
141
            fronts->spaceMyfrontvalue[ACCESSTOFRONT (d - arbitdistances.ins,
 
142
                                                     k + 1)].dptabrow;
 
143
          tmpfront.dptabdirection = MYINSERTIONBIT;
 
144
        }
 
145
      }
 
146
 
 
147
      j = i - k;
 
148
 
 
149
      /* if i = MINUSINFINITYINY or MINUSINFINITYINY + 1 */
 
150
      if (i < 0)
 
151
      {
 
152
        if (tmpfront.dptabdirection == (GtUchar) 0)
 
153
        {
 
154
          alwaysMININFINITYINT = false;
 
155
        }
 
156
        tmpfront.dptabrow = MINUSINFINITYINT;
 
157
      }
 
158
      /* alignment score smaller than T - X */
 
159
      else if ( d_pre > 0 &&
 
160
              big_t.spaceXdropscore != NULL &&
 
161
              (S (i + j, d) < big_t.spaceXdropscore[d_pre] - xdropbelowscore)
 
162
             )
 
163
      {
 
164
        tmpfront.dptabrow = MINUSINFINITYINT;
 
165
      }
 
166
      else if (( k <= -d || k >= d ) || /* not correct boundaries for
 
167
                                           ACCESTOFRONT(d-1,k) */
 
168
               ((fronts->spaceMyfrontvalue[ACCESSTOFRONT(d-1, k)].dptabrow < i)
 
169
                 && (i <= MIN(ulen,vlen + k))))
 
170
      {
 
171
        while (i < ulen && j < vlen)
 
172
        {
 
173
          COMPARESYMBOLSSEP(i,j);
 
174
          i++;
 
175
          j++;
 
176
        }
 
177
 
 
178
        alwaysMININFINITYINT = false;
 
179
        tmpfront.dptabrow = i;
 
180
 
 
181
        /* MAX(bigt_tmp, S(i+j,d)); */
 
182
        if (S (i + j, d) > bigt_tmp)
 
183
        {
 
184
          xdropbest->score = bigt_tmp = S (i + j, d);
 
185
          xdropbest->ivalue = (unsigned int) i;
 
186
          xdropbest->jvalue = (unsigned int) j;
 
187
        }
 
188
      }
 
189
 
 
190
      else
 
191
      {
 
192
        alwaysMININFINITYINT = false;
 
193
        tmpfront.dptabrow =
 
194
          fronts->spaceMyfrontvalue[ACCESSTOFRONT(d - 1, k)].dptabrow;
 
195
      }
 
196
 
 
197
      STOREINARRAYFRONTS (fronts, ACCESSTOFRONT (d, k), Myfrontvalue,
 
198
                            tmpfront);
 
199
 
 
200
      /* delete value for test for */
 
201
      /* INSERTIONBIT/REPLACEMENTBIT/DELETIONBIT above
 
202
      tmpfront.dptabdirection = (GtUchar) 0; */
 
203
 
 
204
    }   /* end for loop */
 
205
 
 
206
    /* if all front values are MINUSINFINITYINT,
 
207
       aligment prematurely finished if allowedMININFINITYINTgenerations
 
208
       exceeded
 
209
       (full front has already ended at d - currentMININFINITYINTgeneration). */
 
210
    if (alwaysMININFINITYINT)
 
211
    {
 
212
      currentMININFINITYINTgeneration++;
 
213
      if (currentMININFINITYINTgeneration > allowedMININFINITYINTgenerations)
 
214
      {
 
215
        d = d - currentMININFINITYINTgeneration;
 
216
        break;
 
217
      }
 
218
    }
 
219
    else
 
220
    {
 
221
      currentMININFINITYINTgeneration = 0;
 
222
      alwaysMININFINITYINT = true;
 
223
    }
 
224
    GT_STOREINARRAY (&big_t, Xdropscore, 10, bigt_tmp);
 
225
 
 
226
    /* fill out of bounds values of MINUSINFINITYINT
 
227
       needed for gt_showmatrix function */
 
228
    for (k = -d; k < lbound - 1; k++)
 
229
    {
 
230
      tmpfront.dptabrow = MINUSINFINITYINT;
 
231
      STOREINARRAYFRONTS (fronts, ACCESSTOFRONT (d, k), Myfrontvalue,
 
232
                          tmpfront);
 
233
    }
 
234
    for (k = ubound + 2; k <= d; k++)
 
235
    {
 
236
      tmpfront.dptabrow = MINUSINFINITYINT;
 
237
      STOREINARRAYFRONTS (fronts, ACCESSTOFRONT (d, k), Myfrontvalue,
 
238
                          tmpfront);
 
239
    }
 
240
 
 
241
    /* alignment finished */
 
242
    if ((-d <= end_k && end_k <= d) &&
 
243
        fronts->spaceMyfrontvalue[ACCESSTOFRONT (d, end_k)].dptabrow == ulen)
 
244
    {
 
245
      break;
 
246
    }
 
247
 
 
248
    /* pruning lower bound
 
249
       lbound may decrease by one or increase/stays the same */
 
250
    for (k = lbound - 1; k <= ubound + 1; k++)
 
251
    {
 
252
      if (fronts->spaceMyfrontvalue[ACCESSTOFRONT (d, k)].dptabrow >
 
253
          MINUSINFINITYINT)
 
254
      {
 
255
        lbound = k;
 
256
        break;
 
257
      }
 
258
    }
 
259
 
 
260
    /* pruning upper bound
 
261
       ubound may increase by one or decrease/stays the same */
 
262
    for (k = ubound + 1; k >= lbound - 1; k--)
 
263
    {
 
264
      if (fronts->spaceMyfrontvalue[ACCESSTOFRONT (d, k)].dptabrow >
 
265
          MINUSINFINITYINT)
 
266
      {
 
267
        ubound = k;
 
268
        break;
 
269
      }
 
270
    }
 
271
 
 
272
    /* handling boundaries lower bound */
 
273
    lboundtmp = lbound;
 
274
    for (k = 0; k >= lbound; k--)
 
275
    {
 
276
      if (fronts->spaceMyfrontvalue[ACCESSTOFRONT (d, k)].dptabrow == vlen + k)
 
277
      {
 
278
        /* debugging 12.11.06, vorher:
 
279
           lboundtmp = k + 2; bei Zhang so */
 
280
        lboundtmp = k;
 
281
        break;
 
282
      }
 
283
    }
 
284
 
 
285
    /* handling boundaries upper bound */
 
286
    uboundtmp = ubound;
 
287
    for (k = 0; k <= ubound; k++)
 
288
    {
 
289
      if (fronts->spaceMyfrontvalue[ACCESSTOFRONT (d, k)].dptabrow == ulen)
 
290
      {
 
291
        /* debugging 12.11.06, vorher:
 
292
           uboundtmp = k - 2; bei Zhang so */
 
293
        uboundtmp = k;
 
294
        break;
 
295
      }
 
296
    }
 
297
 
 
298
    lbound = MAX (lbound, lboundtmp);
 
299
    ubound = MIN (ubound, uboundtmp);
 
300
  }
 
301
 
 
302
  GT_FREEARRAY (&big_t, Xdropscore);
 
303
}