← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/match/tyr-mersplit.c

  • 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/match/tyr-mersplit.c
 
1
/*
 
2
  Copyright (c) 2008 Stefan Kurtz <kurtz@zbh.uni-hamburg.de>
 
3
  Copyright (c) 2008 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
#include <math.h>
 
19
#include <errno.h>
 
20
#include "core/defined-types.h"
 
21
#include "core/divmodmul.h"
 
22
#include "core/fa.h"
 
23
#include "core/intbits.h"
 
24
#include "core/mathsupport.h"
 
25
#include "core/minmax.h"
 
26
#include "core/str_api.h"
 
27
#include "core/types_api.h"
 
28
#include "core/unused_api.h"
 
29
#include "core/xansi_api.h"
 
30
#include "spacedef.h"
 
31
#include "tyr-map.h"
 
32
#include "tyr-mersplit.h"
 
33
 
 
34
#define BUCKETSUFFIX                     ".mbd"
 
35
#define MAXUCHARVALUEWITHBITS(BITNUM)    ((1 << (BITNUM)) - 1)
 
36
#define ISBOUNDDEFINED(UDB,IDX)          GT_ISIBITSET(UDB,IDX)
 
37
#define SETDEFINEDBOUND(UDB,IDX)         GT_SETIBIT(UDB,IDX)
 
38
 
 
39
struct Tyrbckinfo
 
40
{
 
41
  void *mappedmbdfileptr;
 
42
  unsigned int prefixlength;
 
43
  unsigned long numofcodes,
 
44
                *boundisdefined,
 
45
                *bounds;
 
46
  GtUchar remainmask;
 
47
};
 
48
 
 
49
typedef struct
 
50
{
 
51
  const GtUchar *leftmer,
 
52
              *rightmer;
 
53
} Merbounds;
 
54
 
 
55
static unsigned long extractprefixbytecode(unsigned long merbytes,
 
56
                                           unsigned int prefixlength,
 
57
                                           const GtUchar *bytecode)
 
58
{
 
59
  unsigned long idx, code = 0;
 
60
 
 
61
  for (idx=0; idx < MIN((unsigned long) sizeof (unsigned long),merbytes); idx++)
 
62
  {
 
63
    code = (code << 8) | bytecode[idx];
 
64
    if (GT_MULT4(idx+1) == (unsigned long) prefixlength)
 
65
    {
 
66
      break;
 
67
    }
 
68
    if (GT_MULT4(idx+1) > (unsigned long) prefixlength)
 
69
    {
 
70
      code >>= GT_MULT2(GT_MULT4(idx+1) - prefixlength);
 
71
      break;
 
72
    }
 
73
  }
 
74
  return code;
 
75
}
 
76
 
 
77
static GtUchar extractremainingbytes(const GtUchar remainmask,
 
78
                                   unsigned long byteoffset,
 
79
                                   const GtUchar *bytecode)
 
80
{
 
81
  return bytecode[byteoffset] & remainmask;
 
82
}
 
83
 
 
84
static const GtUchar *remainingleftmost(unsigned long merbytes,
 
85
                                      unsigned long byteoffset,
 
86
                                      GtUchar remainmask,
 
87
                                      GtUchar code,
 
88
                                      const GtUchar *leftptr,
 
89
                                      const GtUchar *rightptr)
 
90
{
 
91
  unsigned long len;
 
92
  GtUchar midcode;
 
93
  const GtUchar *midptr;
 
94
 
 
95
  while (leftptr + merbytes < rightptr)
 
96
  {
 
97
    len = (unsigned long) (rightptr-leftptr)/GT_MULT2(merbytes);
 
98
    midptr = leftptr + merbytes * len;
 
99
    midcode = extractremainingbytes(remainmask,byteoffset,midptr);
 
100
    if (code <= midcode)
 
101
    {
 
102
      rightptr = midptr;
 
103
    } else
 
104
    {
 
105
      leftptr = midptr;
 
106
    }
 
107
  }
 
108
  return rightptr;
 
109
}
 
110
 
 
111
static const GtUchar *remainingrightmost(unsigned long merbytes,
 
112
                                       unsigned long byteoffset,
 
113
                                       GtUchar remainmask,
 
114
                                       GtUchar code,
 
115
                                       const GtUchar *leftptr,
 
116
                                       const GtUchar *rightptr)
 
117
{
 
118
  unsigned long len;
 
119
  GtUchar midcode;
 
120
  const GtUchar *midptr;
 
121
 
 
122
  while (leftptr + merbytes < rightptr)
 
123
  {
 
124
    len = (unsigned long) (rightptr-leftptr)/GT_MULT2(merbytes);
 
125
    midptr = leftptr + merbytes * len;
 
126
    midcode = extractremainingbytes(remainmask,byteoffset,midptr);
 
127
    if (code >= midcode)
 
128
    {
 
129
      leftptr = midptr;
 
130
    } else
 
131
    {
 
132
      rightptr = midptr;
 
133
    }
 
134
  }
 
135
  return leftptr;
 
136
}
 
137
 
 
138
static bool remainadvance(Merbounds *merbounds,
 
139
                          unsigned long merbytes,
 
140
                          unsigned long byteoffset,
 
141
                          GtUchar remainmask,
 
142
                          const GtUchar *searchbytecode,
 
143
                          const GtUchar *leftptr,
 
144
                          const GtUchar *rightptr)
 
145
{
 
146
  GtUchar scode, scodeleft, scoderight;
 
147
 
 
148
  scode = extractremainingbytes(remainmask,byteoffset,searchbytecode);
 
149
  scodeleft = extractremainingbytes(remainmask,byteoffset,leftptr);
 
150
  if (scode > scodeleft)
 
151
  {
 
152
    scoderight = extractremainingbytes(remainmask,byteoffset,rightptr);
 
153
    if (scode > scoderight)
 
154
    {
 
155
      return false;
 
156
    }
 
157
    merbounds->leftmer = remainingleftmost(merbytes,
 
158
                                           byteoffset,
 
159
                                           remainmask,
 
160
                                           scode,
 
161
                                           leftptr,
 
162
                                           rightptr);
 
163
  }
 
164
  if (scode < scodeleft)
 
165
  {
 
166
    return false;
 
167
  }
 
168
  if (scode == scodeleft)
 
169
  {
 
170
    merbounds->leftmer = leftptr;
 
171
  }
 
172
  scoderight = extractremainingbytes(remainmask,byteoffset,rightptr);
 
173
  if (scode >= scoderight)
 
174
  {
 
175
    merbounds->rightmer = rightptr;
 
176
  } else
 
177
  {
 
178
    merbounds->rightmer = remainingrightmost(merbytes,
 
179
                                             byteoffset,
 
180
                                             remainmask,
 
181
                                             scode,
 
182
                                             leftptr,
 
183
                                             rightptr);
 
184
  }
 
185
  return true;
 
186
}
 
187
 
 
188
const GtUchar *gt_searchinbuckets(const Tyrindex *tyrindex,
 
189
                             const Tyrbckinfo *tyrbckinfo,
 
190
                             const GtUchar *bytecode)
 
191
{
 
192
  const GtUchar *result;
 
193
  unsigned long prefixcode, leftbound, merbytes;
 
194
 
 
195
  gt_assert(tyrbckinfo != NULL);
 
196
  merbytes = gt_tyrindex_merbytes(tyrindex);
 
197
  prefixcode = extractprefixbytecode(merbytes,
 
198
                                     tyrbckinfo->prefixlength,
 
199
                                     bytecode);
 
200
  leftbound = tyrbckinfo->bounds[prefixcode];
 
201
  if (ISBOUNDDEFINED(tyrbckinfo->boundisdefined,prefixcode))
 
202
  {
 
203
    const GtUchar *mertable;
 
204
    unsigned long rightbound;
 
205
 
 
206
    mertable = gt_tyrindex_mertable(tyrindex);
 
207
    rightbound = tyrbckinfo->bounds[prefixcode+1] - merbytes;
 
208
    if (GT_MOD4(tyrbckinfo->prefixlength) == 0)
 
209
    {
 
210
      result = gt_tyrindex_binmersearch(tyrindex,
 
211
                                     (unsigned long)
 
212
                                     GT_DIV4(tyrbckinfo->prefixlength),
 
213
                                     bytecode,
 
214
                                     mertable + leftbound,
 
215
                                     mertable + rightbound);
 
216
    } else
 
217
    {
 
218
      Merbounds merbounds;
 
219
 
 
220
      merbounds.leftmer = merbounds.rightmer = NULL;
 
221
      if (!remainadvance(&merbounds,
 
222
                         merbytes,
 
223
                         (unsigned long) GT_DIV4(tyrbckinfo->prefixlength),
 
224
                         tyrbckinfo->remainmask,
 
225
                         bytecode,
 
226
                         mertable + leftbound,
 
227
                         mertable + rightbound) ||
 
228
         merbounds.leftmer == NULL ||
 
229
         merbounds.leftmer > merbounds.rightmer)
 
230
      {
 
231
        result = NULL;
 
232
      } else
 
233
      {
 
234
        result = gt_tyrindex_binmersearch(tyrindex,
 
235
                                       1UL + (unsigned long)
 
236
                                             GT_DIV4(tyrbckinfo->prefixlength),
 
237
                                       bytecode,
 
238
                                       merbounds.leftmer,
 
239
                                       merbounds.rightmer);
 
240
      }
 
241
    }
 
242
  } else
 
243
  {
 
244
    result = NULL;
 
245
  }
 
246
  return result;
 
247
}
 
248
 
 
249
static const GtUchar *findrightmostmer(unsigned long merbytes,
 
250
                                     unsigned int prefixlength,
 
251
                                     unsigned long code,
 
252
                                     const GtUchar *leftptr,
 
253
                                     const GtUchar *rightptr)
 
254
{
 
255
  unsigned long len, midcode;
 
256
  const GtUchar *midptr;
 
257
 
 
258
  while (leftptr + merbytes < rightptr)
 
259
  {
 
260
    len = (unsigned long) (rightptr-leftptr)/GT_MULT2(merbytes);
 
261
    midptr = leftptr + merbytes * len;
 
262
    midcode = extractprefixbytecode(merbytes,prefixlength,midptr);
 
263
    if (midcode > code)
 
264
    {
 
265
      rightptr = midptr;
 
266
    } else
 
267
    {
 
268
      leftptr = midptr;
 
269
    }
 
270
  }
 
271
  return leftptr;
 
272
}
 
273
 
 
274
static void splitmerinterval(Tyrbckinfo *tyrbckinfo,
 
275
                             const Tyrindex *tyrindex)
 
276
 
 
277
{
 
278
  const GtUchar *rightbound, *leftptr, *rightptr, *mertable, *lastmer;
 
279
  unsigned long code, leftcode, rightcode, merbytes;
 
280
 
 
281
  mertable = gt_tyrindex_mertable(tyrindex);
 
282
  lastmer = gt_tyrindex_lastmer(tyrindex);
 
283
  merbytes = gt_tyrindex_merbytes(tyrindex);
 
284
  leftptr = mertable;
 
285
  rightptr = lastmer;
 
286
  rightcode = extractprefixbytecode(merbytes,
 
287
                                    tyrbckinfo->prefixlength,
 
288
                                    rightptr);
 
289
  while (true)
 
290
  {
 
291
    leftcode = extractprefixbytecode(merbytes,
 
292
                                     tyrbckinfo->prefixlength,
 
293
                                     leftptr);
 
294
    tyrbckinfo->bounds[leftcode] = (unsigned long) (leftptr - mertable);
 
295
    SETDEFINEDBOUND(tyrbckinfo->boundisdefined,leftcode);
 
296
    if (leftcode == rightcode)
 
297
    {
 
298
      break;
 
299
    }
 
300
    rightbound = findrightmostmer(merbytes,
 
301
                                  tyrbckinfo->prefixlength,
 
302
                                  leftcode,leftptr,rightptr);
 
303
    leftptr = rightbound + merbytes;
 
304
  }
 
305
  tyrbckinfo->bounds[tyrbckinfo->numofcodes]
 
306
    = (unsigned long) (lastmer + merbytes - mertable);
 
307
  SETDEFINEDBOUND(tyrbckinfo->boundisdefined,tyrbckinfo->numofcodes);
 
308
  for (code = tyrbckinfo->numofcodes - 1; /* Nothing */ ; code--)
 
309
  {
 
310
    if (!ISBOUNDDEFINED(tyrbckinfo->boundisdefined,code))
 
311
    {
 
312
      tyrbckinfo->bounds[code] = tyrbckinfo->bounds[code+1];
 
313
    }
 
314
    if (code == 0)
 
315
    {
 
316
      break;
 
317
    }
 
318
  }
 
319
}
 
320
 
 
321
int gt_constructmerbuckets(const char *inputindex,
 
322
                           const Definedunsignedint *callprefixlength,
 
323
                           GtError *err)
 
324
{
 
325
  Tyrindex *tyrindex;
 
326
  Tyrbckinfo tyrbckinfo;
 
327
  FILE *bucketfp = NULL;
 
328
  bool haserr = false;
 
329
 
 
330
  gt_error_check(err);
 
331
  tyrbckinfo.bounds = NULL;
 
332
  tyrbckinfo.boundisdefined = NULL;
 
333
  tyrindex = gt_tyrindex_new(inputindex,err);
 
334
  if (tyrindex == NULL)
 
335
  {
 
336
    haserr = true;
 
337
  }
 
338
  if (!haserr && tyrindex != NULL && !gt_tyrindex_isempty(tyrindex))
 
339
  {
 
340
    if (gt_determinetyrbckpfxlen(&tyrbckinfo.prefixlength,
 
341
                              tyrindex,
 
342
                              callprefixlength,
 
343
                              err) != 0)
 
344
    {
 
345
      haserr = true;
 
346
    }
 
347
  }
 
348
  if (!haserr && tyrindex != NULL && !gt_tyrindex_isempty(tyrindex))
 
349
  {
 
350
    printf("# construct mer buckets for prefixlength %u\n",
 
351
            tyrbckinfo.prefixlength);
 
352
    tyrbckinfo.numofcodes
 
353
      = gt_power_for_small_exponents(gt_tyrindex_alphasize(tyrindex),
 
354
                                     tyrbckinfo.prefixlength);
 
355
    tyrbckinfo.mappedmbdfileptr = NULL;
 
356
    printf("# numofcodes = %lu\n",tyrbckinfo.numofcodes);
 
357
    gt_tyrindex_show(tyrindex);
 
358
    ALLOCASSIGNSPACE(tyrbckinfo.bounds,NULL,unsigned long,
 
359
                     tyrbckinfo.numofcodes+1);
 
360
    GT_INITBITTAB(tyrbckinfo.boundisdefined,tyrbckinfo.numofcodes+1);
 
361
    splitmerinterval(&tyrbckinfo,tyrindex);
 
362
    bucketfp = gt_fa_fopen_with_suffix(inputindex,BUCKETSUFFIX,"wb",err);
 
363
    if (bucketfp == NULL)
 
364
    {
 
365
      haserr = true;
 
366
    }
 
367
  }
 
368
  if (!haserr && tyrindex != NULL && !gt_tyrindex_isempty(tyrindex))
 
369
  {
 
370
    unsigned long pl_long = (unsigned long) tyrbckinfo.prefixlength;
 
371
    gt_assert(bucketfp != NULL);
 
372
    gt_xfwrite(&pl_long, sizeof (pl_long), (size_t) 1, bucketfp);
 
373
  }
 
374
  if (!haserr && tyrindex != NULL && !gt_tyrindex_isempty(tyrindex))
 
375
  {
 
376
    gt_assert(bucketfp != NULL);
 
377
    gt_xfwrite(tyrbckinfo.bounds, sizeof (*tyrbckinfo.bounds),
 
378
               (size_t) (tyrbckinfo.numofcodes+1), bucketfp);
 
379
  }
 
380
  if (!haserr && tyrindex != NULL && !gt_tyrindex_isempty(tyrindex))
 
381
  {
 
382
    gt_assert(bucketfp != NULL);
 
383
    gt_xfwrite(tyrbckinfo.boundisdefined, sizeof (*tyrbckinfo.boundisdefined),
 
384
               GT_NUMOFINTSFORBITS(tyrbckinfo.numofcodes+1), bucketfp);
 
385
  }
 
386
  gt_fa_xfclose(bucketfp);
 
387
  if (tyrindex != NULL)
 
388
  {
 
389
    gt_tyrindex_delete(&tyrindex);
 
390
  }
 
391
  FREESPACE(tyrbckinfo.bounds);
 
392
  gt_free(tyrbckinfo.boundisdefined);
 
393
  return haserr ? -1 : 0;
 
394
}
 
395
 
 
396
Tyrbckinfo *gt_tyrbckinfo_new(const char *tyrindexname,unsigned int alphasize,
 
397
                              GtError *err)
 
398
{
 
399
  size_t numofbytes;
 
400
  Tyrbckinfo *tyrbckinfo;
 
401
  bool haserr = false;
 
402
 
 
403
  ALLOCASSIGNSPACE(tyrbckinfo,NULL,Tyrbckinfo,1);
 
404
  tyrbckinfo->mappedmbdfileptr = gt_fa_mmap_read_with_suffix(tyrindexname,
 
405
                                                          BUCKETSUFFIX,
 
406
                                                          &numofbytes,err);
 
407
  if (tyrbckinfo->mappedmbdfileptr == NULL)
 
408
  {
 
409
    haserr = true;
 
410
  }
 
411
  if (!haserr)
 
412
  {
 
413
    unsigned long pl_long;
 
414
 
 
415
    gt_assert(tyrbckinfo->mappedmbdfileptr != NULL);
 
416
    pl_long = *((unsigned long *) tyrbckinfo->mappedmbdfileptr);
 
417
    tyrbckinfo->prefixlength = (unsigned int) pl_long;
 
418
    tyrbckinfo->numofcodes
 
419
      = gt_power_for_small_exponents(alphasize,tyrbckinfo->prefixlength);
 
420
    /*check if numofbytes == expected size*/
 
421
    gt_assert(numofbytes == sizeof (unsigned long) *
 
422
                            (1UL +
 
423
                             (tyrbckinfo->numofcodes+1) +
 
424
                             GT_NUMOFINTSFORBITS(tyrbckinfo->numofcodes + 1)));
 
425
    tyrbckinfo->bounds = ((unsigned long *) tyrbckinfo->mappedmbdfileptr) + 1;
 
426
    tyrbckinfo->boundisdefined
 
427
      = tyrbckinfo->bounds + tyrbckinfo->numofcodes + 1;
 
428
    if (tyrbckinfo->prefixlength > 0 && GT_MOD4(tyrbckinfo->prefixlength) > 0)
 
429
    {
 
430
      tyrbckinfo->remainmask
 
431
        = (GtUchar) MAXUCHARVALUEWITHBITS(GT_MULT2(
 
432
                                       4U - GT_MOD4(tyrbckinfo->prefixlength)));
 
433
    }
 
434
  }
 
435
  if (haserr)
 
436
  {
 
437
    FREESPACE(tyrbckinfo);
 
438
    return NULL;
 
439
  }
 
440
  return tyrbckinfo;
 
441
}
 
442
 
 
443
void gt_tyrbckinfo_delete(Tyrbckinfo **tyrbckinfoptr)
 
444
{
 
445
  Tyrbckinfo *tyrbckinfo = *tyrbckinfoptr;
 
446
 
 
447
  gt_fa_xmunmap(tyrbckinfo->mappedmbdfileptr);
 
448
  tyrbckinfo->mappedmbdfileptr = NULL;
 
449
  FREESPACE(tyrbckinfo);
 
450
  *tyrbckinfoptr = NULL;
 
451
}