source: branches/lib/SL/FILTER/AP_filter.cxx

Last change on this file was 17877, checked in by westram, 7 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.4 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : AP_filter.cxx                                     //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include "AP_filter.hxx"
12#include <arbdb.h>
13
14// ------------------
15//      AP_filter
16
17void AP_filter::init(size_t size) {
18    filter_mask        = new bool[size];
19    filter_len         = size;
20    real_len           = 0;
21    update             = AP_timer();
22    simplify_type      = AWT_FILTER_SIMPLIFY_NOT_INITIALIZED;
23    simplify[0]        = 0;   // silence cppcheck-warning
24    bootstrap          = NULp;
25    filterpos_2_seqpos = NULp;
26
27#if defined(ASSERTION_USED)
28    checked_for_validity = false;
29#endif
30}
31
32void AP_filter::make_permeable(size_t size) {
33    init(size);
34    real_len = filter_len;
35    for (size_t i = 0; i < size; i++) filter_mask[i] = true;
36}
37
38void AP_filter::init_from_string(const char *ifilter, const char *zerobases, size_t size) {
39    init(size);
40
41    bool   char2mask[256];
42    size_t i;
43
44    for (i = 0; i<256; ++i) char2mask[i] = true;
45    if (zerobases) {
46        for (i = 0; zerobases[i]; ++i) char2mask[safeCharIndex(zerobases[i])] = false;
47    }
48    else {
49        char2mask['0'] = false;
50    }
51
52    real_len = 0;
53    for (i = 0; i < size && ifilter[i]; ++i) {
54        real_len += int(filter_mask[i] = char2mask[safeCharIndex(ifilter[i])]);
55    }
56    for (; i < size; i++) {
57        filter_mask[i] = true;
58        real_len++;
59    }
60}
61
62
63AP_filter::AP_filter(size_t size) {
64    make_permeable(size);
65}
66
67AP_filter::AP_filter(const AP_filter& other)
68    : filter_mask(new bool[other.filter_len]),
69      filter_len(other.filter_len),
70      real_len(other.real_len),
71      update(other.update),
72      simplify_type(other.simplify_type),
73      bootstrap(NULp),
74      filterpos_2_seqpos(NULp)
75{
76    memcpy(filter_mask, other.filter_mask, filter_len*sizeof(*filter_mask));
77    memcpy(simplify, other.simplify, sizeof(simplify)*sizeof(*simplify));
78    if (other.bootstrap) {
79        bootstrap = new size_t[real_len];
80        memcpy(bootstrap, other.bootstrap, real_len*sizeof(*bootstrap));
81    }
82    if (other.filterpos_2_seqpos) {
83        filterpos_2_seqpos = new size_t[real_len];
84        memcpy(filterpos_2_seqpos, other.filterpos_2_seqpos, real_len*sizeof(*filterpos_2_seqpos));
85    }
86#if defined(ASSERTION_USED)
87    checked_for_validity = other.checked_for_validity;
88#endif
89}
90
91
92
93AP_filter::AP_filter(const char *ifilter, const char *zerobases, size_t size) {
94    if (!ifilter || !*ifilter) {
95        make_permeable(size);
96    }
97    else {
98        init_from_string(ifilter, zerobases, size);
99    }
100}
101
102AP_filter::AP_filter(AF_Not, const AP_filter& other) {
103    size_t      size  = other.get_length();
104    const bool *omask = other.filter_mask;
105
106    init(size);
107    for (size_t i = 0; i < size; i++) {
108        real_len += (filter_mask[i] = !omask[i]);
109    }
110}
111
112AP_filter::AP_filter(const AP_filter& f1, AF_Combine comb, const AP_filter& f2) {
113    size_t size = f1.get_length();
114    af_assert(size == f2.get_length());
115
116    init(size);
117
118    const bool *m1 = f1.filter_mask;
119    const bool *m2 = f2.filter_mask;
120
121    switch (comb) {
122        case AND:
123            for (size_t i = 0; i<size; ++i) {
124                real_len += (filter_mask[i] = (m1[i] && m2[i]));
125            }
126            break;
127        case OR:
128            for (size_t i = 0; i<size; ++i) {
129                real_len += (filter_mask[i] = (m1[i] || m2[i]));
130            }
131            break;
132        case XOR:
133            for (size_t i = 0; i<size; ++i) {
134                real_len += (filter_mask[i] = (m1[i] ^ m2[i]));
135            }
136            break;
137    }
138}
139
140AP_filter::~AP_filter() {
141    delete [] bootstrap;
142    delete [] filter_mask;
143    delete [] filterpos_2_seqpos;
144}
145
146char *AP_filter::to_string() const {
147    af_assert(checked_for_validity);
148
149    char *data = ARB_alloc<char>(filter_len+1);
150
151    for (size_t i=0; i<filter_len; ++i) {
152        data[i] = "01"[filter_mask[i]];
153    }
154    data[filter_len] = 0;
155
156    return data;
157}
158
159
160void AP_filter::enable_simplify(AWT_FILTER_SIMPLIFY type) {
161    if (type != simplify_type) {
162        int i;
163        for (i=0; i<32; i++) {
164            simplify[i] = '.';
165        }
166        for (; i<256; i++) { // LOOP_VECTORIZED // tested down to gcc 5.5.0 (may fail on older gcc versions)
167            simplify[i] = i;
168        }
169        switch (type) {
170            case AWT_FILTER_SIMPLIFY_DNA:
171                simplify[(unsigned char)'g'] = 'a';
172                simplify[(unsigned char)'G'] = 'A';
173                simplify[(unsigned char)'u'] = 'c';
174                simplify[(unsigned char)'t'] = 'c';
175                simplify[(unsigned char)'U'] = 'C';
176                simplify[(unsigned char)'T'] = 'C';
177                break;
178            case AWT_FILTER_SIMPLIFY_PROTEIN:
179                af_assert(0);                           // not implemented or impossible!?
180                break;
181            case AWT_FILTER_SIMPLIFY_NONE:
182                break;
183            default:
184                af_assert(0);
185                break;
186        }
187
188        simplify_type = type;
189    }
190}
191
192void AP_filter::calc_filterpos_2_seqpos() {
193    af_assert(checked_for_validity);
194    af_assert(real_len>0);
195
196    delete [] filterpos_2_seqpos;
197    filterpos_2_seqpos = new size_t[real_len];
198    size_t i, j;
199    for (i=j=0; i<filter_len; ++i) {
200        if (filter_mask[i]) {
201            filterpos_2_seqpos[j++] = i;
202        }
203    }
204}
205
206void AP_filter::enable_bootstrap() {
207    af_assert(checked_for_validity);
208    af_assert(real_len>0);
209
210    delete [] bootstrap;
211    bootstrap = new size_t[real_len];
212
213    af_assert(filter_len < RAND_MAX);
214
215    for (size_t i = 0; i<real_len; ++i) {
216        int r = GB_random(real_len);
217        af_assert(r >= 0);     // otherwise overflow in random number generator
218        bootstrap[i] = r;
219    }
220}
221
222char *AP_filter::blowup_string(const char *filtered_string, char fillChar) const {
223    /*! blow up 'filtered_string' to unfiltered length
224     * by inserting 'fillChar' at filtered positions
225     */
226    af_assert(checked_for_validity);
227
228    char   *blownup = ARB_alloc<char>(filter_len+1);
229    size_t  f       = 0;
230
231    for (size_t i = 0; i<filter_len; ++i) {
232        blownup[i] = use_position(i) ? filtered_string[f++] : fillChar;
233    }
234    blownup[filter_len] = 0;
235
236    return blownup;
237}
238
239char *AP_filter::filter_string(const char *fulllen_string) const {
240    /*! filter given 'fulllen_string'
241     */
242
243    af_assert(checked_for_validity);
244
245    char   *filtered = ARB_alloc<char>(real_len+1);
246    size_t  f        = 0;
247
248    get_filterpos_2_seqpos(); // create if missing
249    for (size_t i = 0; i<real_len; ++i) {
250        size_t p      = filterpos_2_seqpos[i];
251        filtered[f++] = fulllen_string[p];
252    }
253    filtered[f] = 0;
254
255    return filtered;
256}
257
258// -------------------
259//      AP_weights
260
261AP_weights::AP_weights(const AP_filter *fil)
262    : len(fil->get_filtered_length()),
263      weights(NULp)
264{
265}
266
267AP_weights::AP_weights(const GB_UINT4 *w, size_t wlen, const AP_filter *fil)
268    : len(fil->get_filtered_length()),
269      weights(NULp)
270{
271    ARB_alloc_aligned(weights, len);
272
273    af_assert(wlen == fil->get_length());
274
275    size_t i, j;
276    for (j=i=0; j<wlen; ++j) {
277        if (fil->use_position(j)) {
278            weights[i++] = w[j];
279        }
280    }
281    af_assert(j <= fil->get_length());
282    af_assert(i == fil->get_filtered_length());
283}
284
285AP_weights::AP_weights(const AP_weights& other)
286    : len(other.len),
287      weights(NULp)
288{
289    if (other.weights) {
290        ARB_alloc_aligned(weights, len);
291        memcpy(weights, other.weights, len*sizeof(*weights));
292    }
293}
294
295AP_weights::~AP_weights() {
296    free(weights);
297}
298
299long AP_timer() {
300    static long time = 0;
301    return ++time;
302}
303
304
305// --------------------------------------------------------------------------------
306
307#ifdef UNIT_TESTS
308#ifndef TEST_UNIT_H
309#include <test_unit.h>
310#endif
311
312#define TEST_EXPECT_EQUAL_FILTERS(f1,f2) do{            \
313        TEST_EXPECT_NO_ERROR((f1).is_invalid());        \
314        TEST_EXPECT_NO_ERROR((f2).is_invalid());        \
315        char *m1 = (f1).to_string();                    \
316        char *m2 = (f2).to_string();                    \
317        TEST_EXPECT_EQUAL(m1, m2);                      \
318        free(m2);                                       \
319        free(m1);                                       \
320    }while(0)
321
322void TEST_filter() {
323    const int   LEN            = 20;
324    const int   MASK_BITCOUNT  = 9;
325    const char *mask           = "01100001110000110011";
326    const char *mask_inv       = "10011110001111001100";
327    const char *mask_some      = "00100101100101011001";
328    const char *seq            = "MSKTAYTKVLFDRGSALDGK";
329    const char *seq_masked     =  "SK"  "KVL"  "SA""GK";
330    const char *blow_mask      = "_SK____KVL____SA__GK";
331    const char *seq_masked_inv = "M""TAYT" "FDRG""LD";
332    const char *blow_mask_inv  = "M__TAYT___FDRG__LD__";
333
334    AP_filter f1(LEN);
335    AP_filter f2(mask, "0", LEN);
336    AP_filter n2(mask, "1", LEN);
337    AP_filter f3(mask_inv, "0", LEN);
338    AP_filter n3(mask_inv, "1", LEN);
339
340    TEST_EXPECT_EQUAL(f1.get_length(), LEN);
341    TEST_EXPECT_EQUAL(f2.get_length(), LEN);
342    TEST_EXPECT_EQUAL(f3.get_length(), LEN);
343    TEST_EXPECT_EQUAL(n2.get_length(), LEN);
344    TEST_EXPECT_EQUAL(n3.get_length(), LEN);
345
346    TEST_EXPECT_EQUAL(f1.get_filtered_length(), LEN);
347    TEST_EXPECT_EQUAL(f2.get_filtered_length(), MASK_BITCOUNT);
348    TEST_EXPECT_EQUAL(f3.get_filtered_length(), LEN-MASK_BITCOUNT);
349    TEST_EXPECT_EQUAL(n2.get_filtered_length(), LEN-MASK_BITCOUNT);
350    TEST_EXPECT_EQUAL(n3.get_filtered_length(), MASK_BITCOUNT);
351
352    TEST_EXPECT_NO_ERROR(f1.is_invalid());
353    TEST_EXPECT_NO_ERROR(f2.is_invalid());
354    TEST_EXPECT_NO_ERROR(f3.is_invalid());
355    TEST_EXPECT_NO_ERROR(n2.is_invalid());
356    TEST_EXPECT_NO_ERROR(n3.is_invalid());
357
358    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(f1.to_string(), "11111111111111111111");
359    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(f2.to_string(), mask);
360    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(f3.to_string(), mask_inv);
361    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(n2.to_string(), mask_inv);
362    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(n3.to_string(), mask);
363
364    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(f1.filter_string(seq), seq);
365    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(f2.filter_string(seq), seq_masked);
366    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(f3.filter_string(seq), seq_masked_inv);
367    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(n2.filter_string(seq), seq_masked_inv);
368    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(n3.filter_string(seq), seq_masked);
369
370    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(f1.blowup_string(seq,            '_'), seq);
371    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(f2.blowup_string(seq_masked,     '_'), blow_mask);
372    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(f3.blowup_string(seq_masked_inv, '_'), blow_mask_inv);
373
374    // test inverting filters:
375    AP_filter i2(NOT, f2);
376    AP_filter i3(NOT, f3);
377
378    TEST_EXPECT_EQUAL_FILTERS(i2, n2);
379    TEST_EXPECT_EQUAL_FILTERS(i3, n3);
380
381    // test filter combination (AND + OR):
382    AP_filter s2(mask_some, "0", LEN);
383    AP_filter s3(NOT, s2);
384
385    AP_filter as23(s2, AND, s3);
386    AP_filter os23(s2, OR,  s3);
387
388    TEST_EXPECT_ERROR_CONTAINS(as23.is_invalid(), "Sequence completely filtered out (no columns left)");
389
390    TEST_EXPECT_EQUAL(as23.get_filtered_length(), 0);
391    TEST_EXPECT_EQUAL(os23.get_filtered_length(), LEN);
392
393    AP_filter fs22(f2, AND, s2);
394    AP_filter fs23(f2, AND, s3);
395    AP_filter o2223(fs22, OR, fs23);
396
397    TEST_EXPECT_EQUAL_FILTERS(o2223, f2);
398
399    AP_filter x(fs22, XOR, fs23);
400    AP_filter xa1(AP_filter(fs22, AND, AP_filter(NOT, fs23)), OR,  AP_filter(AP_filter(NOT, fs22), AND, fs23)); // = (a&&!b) || (!a&&b)
401    AP_filter xa2(AP_filter(fs22, OR,  fs23),                 AND, AP_filter(NOT, AP_filter(fs22, AND, fs23))); // = (a||b) && !(a&&b)
402
403    TEST_EXPECT_EQUAL_FILTERS(x, xa1);
404    TEST_EXPECT_EQUAL_FILTERS(x, xa2);
405
406    TEST_EXPECT_EQUAL_STRINGCOPY__NOERROREXPORTED(x.to_string(), mask);
407}
408
409#endif // UNIT_TESTS
410
411// --------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.