← Back to branch summary

~ubuntu-branches/debian/jessie/arb/jessie

« back to all changes in this revision

Viewing changes to PROBE_SET/ps_detect_weak_differences.cxx

  • Committer: Package Import Robot
  • Author(s): Elmar Pruesse, Andreas Tille, Elmar Pruesse
  • Date: 2014-09-02 15:15:06 UTC
  • mfrom: (1.1.6)
  • Revision ID: package-import@ubuntu.com-20140902151506-jihq58b3iz342wif
Tags: 6.0.2-1
http://bugs.debian.org/741890
http://bugs.debian.org/757497
[ Andreas Tille ]
* New upstream version
  Closes: #741890
* debian/upstream -> debian/upstream/metadata
* debian/control:
   - Build-Depends: added libglib2.0-dev
   - Depends: added mafft, mrbayes
* debian/rules
   - Add explicite --remove-section=.comment option to manual strip call
* cme fix dpkg-control
* arb-common.dirs: Do not create unneeded lintian dir
* Add turkish debconf translation (thanks for the patch to Mert Dirik
  <mertdirik@gmail.com>)
  Closes: #757497

[ Elmar Pruesse ]
* patches removed:
   - 10_config.makefiles.patch,
     80_no_GL.patch
       removed in favor of creating file from config.makefile.template via 
       sed in debian/control
   - 20_Makefile_main.patch
       merged upstream
   - 21_Makefiles.patch
       no longer needed
   - 30_tmpfile_CVE-2008-5378.patch: 
       merged upstream
   - 50_fix_gcc-4.8.patch:
       merged upstream
   - 40_add_libGLU.patch:
       libGLU not needed for arb_ntree)
   - 60_use_debian_packaged_raxml.patch:
       merged upstream
   - 70_hardening.patch
       merged upstream
   - 72_add_math_lib_to_linker.patch
       does not appear to be needed
* patches added:
   - 10_upstream_r12793__show_db_load_progress:
       backported patch showing progress while ARB is loading a database
       (needed as indicator/splash screen while ARB is launching)
   - 20_upstream_r12794__socket_permissions:
       backported security fix
   - 30_upstream_r12814__desktop_keywords:
       backported add keywords to desktop (fixes lintian warning)
   - 40_upstream_r12815__lintian_spelling:
       backported fix for lintian reported spelling errors
   - 50_private_nameservers
       change configuration to put nameservers into users home dirs
       (avoids need for shared writeable directory)
   - 60_use_debian_phyml
       use phyml from debian package for both interfaces in ARB
* debian/rules:
   - create config.makefile from override_dh_configure target
   - use "make tarfile" in override_dh_install
   - remove extra cleaning not needed for ARB 6
   - use "dh_install --list-missing" to avoid missing files
   - added override_dh_fixperms target
* debian/control:
   - added libarb-dev package
   - Depends: added phyml, xdg-utils
   - Suggests: removed phyml
   - fix lintian duplicate-short-description (new descriptions)
* debian/*.install:
   - "unrolled" confusing globbing to select files
   - pick files from debian/tmp
   - moved all config files to /etc/arb
* debian/arb-common.templates: updated
* scripts:
   - removed arb-add-pt-server
   - launch-wrapper: 
     - only add demo.arb to newly created $ARBUSERDATA
     - pass commandline arguments through bin/arb wrapper
   - preinst: removing old PT server index files on upgrade from 5.5*
   - postinst: set setgid on shared PT dir
* rewrote arb.1 manfile
* added file icon for ARB databases
* using upstream arb_tcp.dat

Show diffs side-by-side

added added

removed removed

Lines of Context:
PROBE_SET/ps_detect_weak_differences.cxx
 
1
// =============================================================== //
 
2
//                                                                 //
 
3
//   File      : ps_detect_weak_differences.cxx                    //
 
4
//   Purpose   :                                                   //
 
5
//                                                                 //
 
6
//   Coded by Wolfram Foerster in October 2002                     //
 
7
//   Institute of Microbiology (Technical University Munich)       //
 
8
//   http://www.arb-home.de/                                       //
 
9
//                                                                 //
 
10
// =============================================================== //
 
11
 
1
12
#include "ps_database.hxx"
2
13
#include "ps_bitmap.hxx"
3
14
#include "ps_tools.hxx"
4
15
 
5
 
#include <cstdio>
6
 
#include <cstdlib>
7
 
#include <cstring>
8
16
#include <climits>
9
 
 
10
17
#include <sys/times.h>
11
18
 
12
19
// common globals
13
 
SpeciesID          __MAX_ID;
14
 
SpeciesID          __MIN_ID;
15
 
PS_BitMap_Fast    *__MAP;
 
20
static SpeciesID          __MAX_ID;
 
21
static SpeciesID          __MIN_ID;
 
22
static PS_BitMap_Fast    *__MAP;
16
23
// globals for PS_detect_weak_differences
17
 
IDVector          *__PATH;
18
 
IDVector          *__INVERSE_PATH;
19
 
unsigned long int  __COUNT_SET_OPS  = 0;
20
 
unsigned long int  __COUNT_SET_OPS2 = 0;
21
 
char              *__NODES_LEFT;
 
24
static IDVector          *__PATH;
 
25
static IDVector          *__INVERSE_PATH;
 
26
static unsigned long int  __COUNT_SET_OPS  = 0;
 
27
static unsigned long int  __COUNT_SET_OPS2 = 0;
 
28
static char              *__NODES_LEFT;
22
29
// globals for PS_print_and_evaluate_map
23
 
IDSet             *__PATHSET;
24
 
IDID2IDSetMap     *__PAIR2PATH;
25
 
SpeciesID          __ONEMATCH_MIN_ID;
26
 
SpeciesID          __ONEMATCH_MAX_ID;
27
 
 
28
 
void PS_print_path() {
29
 
    printf( "__PATH %3zu :",__PATH->size() );
30
 
    int c = 1;
31
 
    for (IDVectorCIter i = __PATH->begin(); i != __PATH->end(); ++i,++c) {
32
 
        if (c % 20 == 0) printf( "\n" );
33
 
        printf( " %3i", *i );
34
 
    }
35
 
    printf( "\n" );
36
 
}
37
 
 
38
 
void PS_print_inverse_path() {
39
 
    printf( "__INVERSE_PATH %3zu :",__INVERSE_PATH->size() );
40
 
    int c = 1;
41
 
    for (IDVectorCIter i = __INVERSE_PATH->begin(); i != __INVERSE_PATH->end(); ++i,++c) {
42
 
        if (c % 20 == 0) printf( "\n" );
43
 
        printf( " %3i", *i );
44
 
    }
45
 
    printf( "\n" );
46
 
}
47
 
 
48
 
 
49
 
void PS_detect_weak_differences_stepdown( const PS_NodePtr _ps_node,
50
 
                                          const SpeciesID  _parent_ID,
51
 
                                          const long       _depth ) {
 
30
static IDSet             *__PATHSET;
 
31
static IDID2IDSetMap     *__PAIR2PATH;
 
32
static SpeciesID          __ONEMATCH_MIN_ID;
 
33
static SpeciesID          __ONEMATCH_MAX_ID;
 
34
 
 
35
#if defined(DEBUG)
 
36
void PS_debug_print_path() {
 
37
    printf("__PATH %3zu :", __PATH->size());
 
38
    int c = 1;
 
39
    for (IDVectorCIter i = __PATH->begin(); i != __PATH->end(); ++i, ++c) {
 
40
        if (c % 20 == 0) printf("\n");
 
41
        printf(" %3i", *i);
 
42
    }
 
43
    printf("\n");
 
44
}
 
45
 
 
46
void PS_debug_print_inverse_path() {
 
47
    printf("__INVERSE_PATH %3zu :", __INVERSE_PATH->size());
 
48
    int c = 1;
 
49
    for (IDVectorCIter i = __INVERSE_PATH->begin(); i != __INVERSE_PATH->end(); ++i, ++c) {
 
50
        if (c % 20 == 0) printf("\n");
 
51
        printf(" %3i", *i);
 
52
    }
 
53
    printf("\n");
 
54
}
 
55
#endif
 
56
 
 
57
static void PS_detect_weak_differences_stepdown(const PS_NodePtr& _ps_node,
 
58
                                                const SpeciesID  _parent_ID,
 
59
                                                const long       _depth) {
52
60
    //  Recursively walk through tree and make a bool-matrix of SpeciesID's
53
61
    //  where true means that the 2 species can be distinguished by a probe.
54
62
    //
58
66
 
59
67
    SpeciesID id = _ps_node->getNum();
60
68
    if (_depth < 60) {
61
 
        printf( "%s", __NODES_LEFT );
62
 
        for (int i = 0; i < 60; ++i) printf( "\b" );
63
 
        fflush( stdout );
 
69
        printf("%s", __NODES_LEFT);
 
70
        for (int i = 0; i < 60; ++i) printf("\b");
 
71
        fflush(stdout);
64
72
    }
65
73
    //
66
74
    // append IDs to paths
67
75
    //
68
 
    __PATH->push_back( id );                                                                    // append id to path
 
76
    __PATH->push_back(id);                                                                      // append id to path
69
77
    for (SpeciesID i = (_parent_ID < __MIN_ID) ? __MIN_ID : _parent_ID+1; i < id; ++i) {        // append parent_id+1 .. id-1 to inverse path
70
 
        //printf( "%i ",i );
71
 
        __INVERSE_PATH->push_back( i );
 
78
        __INVERSE_PATH->push_back(i);
72
79
    }
73
80
 
74
81
    //
76
83
    //
77
84
    if (_ps_node->hasProbes()) {
78
85
        if (_ps_node->hasPositiveProbes() && _ps_node->hasInverseProbes()) {
79
 
//            PS_print_path();
80
 
//            PS_print_inverse_path();
81
86
            unsigned long int set_ops = 2*__PATH->size()*(__MAX_ID-id-1+__INVERSE_PATH->size());
82
87
            if (ULONG_MAX - __COUNT_SET_OPS < set_ops) {
83
88
                set_ops = set_ops - (ULONG_MAX - __COUNT_SET_OPS);
95
100
                // inverse path loop (explicit)
96
101
                for (IDVectorCIter it_inverse_path = __INVERSE_PATH->begin();
97
102
                     it_inverse_path != __INVERSE_PATH->end();
98
 
                     ++it_inverse_path ) {
 
103
                     ++it_inverse_path) {
99
104
                    inverse_path_ID = *it_inverse_path;
100
105
 
101
 
                    __MAP->setTrue( path_ID, inverse_path_ID );
102
 
                    __MAP->setTrue( inverse_path_ID, path_ID );
 
106
                    __MAP->setTrue(path_ID, inverse_path_ID);
 
107
                    __MAP->setTrue(inverse_path_ID, path_ID);
103
108
                }
104
109
 
105
110
                // inverse path loop (implicit)
106
111
                for (inverse_path_ID = id+1; inverse_path_ID < __MAX_ID; ++inverse_path_ID) { // skip to id ABOVE current node id
107
 
                    __MAP->setTrue( path_ID, inverse_path_ID );
108
 
                    __MAP->setTrue( inverse_path_ID, path_ID );
 
112
                    __MAP->setTrue(path_ID, inverse_path_ID);
 
113
                    __MAP->setTrue(inverse_path_ID, path_ID);
109
114
                }
110
115
            }
111
 
        } else {
112
 
//            PS_print_path();
113
 
//            PS_print_inverse_path();
 
116
        }
 
117
        else {
114
118
            unsigned long int set_ops = __PATH->size()*(__MAX_ID-id-1+__INVERSE_PATH->size());
115
119
            if (ULONG_MAX - __COUNT_SET_OPS < set_ops) {
116
120
                set_ops = set_ops - (ULONG_MAX - __COUNT_SET_OPS);
130
134
                // inverse path loop (explicit)
131
135
                for (IDVectorCIter it_inverse_path = __INVERSE_PATH->begin();
132
136
                     it_inverse_path != __INVERSE_PATH->end();
133
 
                     ++it_inverse_path ) {
 
137
                     ++it_inverse_path) {
134
138
                    inverse_path_ID = *it_inverse_path;
135
139
                    smaller_ID = (path_ID < inverse_path_ID) ? path_ID : inverse_path_ID;
136
140
                    bigger_ID  = (path_ID > inverse_path_ID) ? path_ID : inverse_path_ID;
137
141
 
138
 
                    if (__MAP->get( smaller_ID, bigger_ID )) {
139
 
                        __MAP->setTrue( bigger_ID, smaller_ID );
140
 
                    } else {
141
 
                        __MAP->setTrue( smaller_ID, bigger_ID );
 
142
                    if (__MAP->get(smaller_ID, bigger_ID)) {
 
143
                        __MAP->setTrue(bigger_ID, smaller_ID);
 
144
                    }
 
145
                    else {
 
146
                        __MAP->setTrue(smaller_ID, bigger_ID);
142
147
                    }
143
148
                }
144
149
 
147
152
                    smaller_ID = (path_ID < inverse_path_ID) ? path_ID : inverse_path_ID;
148
153
                    bigger_ID  = (path_ID > inverse_path_ID) ? path_ID : inverse_path_ID;
149
154
 
150
 
                    if (__MAP->get( smaller_ID, bigger_ID )) {
151
 
                        __MAP->setTrue( bigger_ID, smaller_ID );
152
 
                    } else {
153
 
                        __MAP->setTrue( smaller_ID, bigger_ID );
 
155
                    if (__MAP->get(smaller_ID, bigger_ID)) {
 
156
                        __MAP->setTrue(bigger_ID, smaller_ID);
 
157
                    }
 
158
                    else {
 
159
                        __MAP->setTrue(smaller_ID, bigger_ID);
154
160
                    }
155
161
                }
156
162
            }
160
166
    // step down the children
161
167
    //
162
168
    int c = _ps_node->countChildren()-1;
163
 
    for ( PS_NodeMapConstIterator i = _ps_node->getChildrenBegin();
 
169
    for (PS_NodeMapConstIterator i = _ps_node->getChildrenBegin();
164
170
          i != _ps_node->getChildrenEnd();
165
 
          ++i,--c ) {
 
171
          ++i, --c) {
166
172
        if (_depth < 60) {
167
173
            if (c < 10) {
168
 
                __NODES_LEFT[ _depth ] = '0'+c;
169
 
            } else {
170
 
                __NODES_LEFT[ _depth ] = '+';
 
174
                __NODES_LEFT[_depth] = '0'+c;
 
175
            }
 
176
            else {
 
177
                __NODES_LEFT[_depth] = '+';
171
178
            }
172
179
        }
173
 
        PS_detect_weak_differences_stepdown( i->second, id, _depth+1 );
 
180
        PS_detect_weak_differences_stepdown(i->second, id, _depth+1);
174
181
    }
175
 
    if (_depth < 60) __NODES_LEFT[ _depth ] = ' ';
 
182
    if (_depth < 60) __NODES_LEFT[_depth] = ' ';
176
183
 
177
184
    //
178
185
    // remove IDs from paths
183
190
    }
184
191
}
185
192
 
186
 
void PS_detect_weak_differences( const PS_NodePtr _root_node ) {
 
193
static void PS_detect_weak_differences(const PS_NodePtr& _root_node) {
187
194
    //
188
195
    // make bitmap
189
196
    //
192
199
 
193
200
    int c = 0;
194
201
    struct tms before;
195
 
    times( &before );
 
202
    times(&before);
196
203
    struct tms before_first_level_node;
197
 
    for (PS_NodeMapConstIterator i = _root_node->getChildrenBegin(); i != _root_node->getChildrenEnd(); ++i,++c ) {
 
204
    for (PS_NodeMapConstIterator i = _root_node->getChildrenBegin(); i != _root_node->getChildrenEnd(); ++i, ++c) {
198
205
        if (_root_node->countChildren()-c-1 < 10) {
199
206
            __NODES_LEFT[0] = '0'+_root_node->countChildren()-c-1;
200
 
        } else {
 
207
        }
 
208
        else {
201
209
            __NODES_LEFT[0] = '+';
202
210
        }
203
211
        if ((c < 50) || (c % 100 == 0)) {
204
 
            times( &before_first_level_node );
205
 
            printf( "PS_detect_weak_differences_stepdown( %i ) : %i. of %zu  ", i->first, c+1, _root_node->countChildren() ); fflush( stdout );
 
212
            times(&before_first_level_node);
 
213
            printf("PS_detect_weak_differences_stepdown( %i ) : %i. of %zu  ", i->first, c+1, _root_node->countChildren()); fflush(stdout);
206
214
        }
207
 
        PS_detect_weak_differences_stepdown( i->second, -1, 1 );
 
215
        PS_detect_weak_differences_stepdown(i->second, -1, 1);
208
216
        if ((c < 50) || (c % 100 == 0)) {
209
 
            PS_print_time_diff( &before_first_level_node, "this node ", "  " );
210
 
            PS_print_time_diff( &before, "total ", "\n" );
 
217
            PS_print_time_diff(&before_first_level_node, "this node ", "  ");
 
218
            PS_print_time_diff(&before, "total ", "\n");
211
219
        }
212
220
    }
213
 
    printf( "%lu * %lu + %lu set operations performed\n", __COUNT_SET_OPS2, ULONG_MAX, __COUNT_SET_OPS ); fflush( stdout );
 
221
    printf("%lu * %lu + %lu set operations performed\n", __COUNT_SET_OPS2, ULONG_MAX, __COUNT_SET_OPS); fflush(stdout);
214
222
 
215
223
    delete __PATH;
216
224
    delete __INVERSE_PATH;
217
225
}
218
226
 
219
 
typedef map<ID2IDPair,PS_NodePtr>    IDID2NodeMap;
 
227
typedef map<ID2IDPair, PS_NodePtr>    IDID2NodeMap;
220
228
typedef IDID2NodeMap::iterator       IDID2NodeMapIter;
221
229
typedef IDID2NodeMap::const_iterator IDID2NodeMapCIter;
222
230
 
223
 
void PS_find_probes_for_pairs( const PS_NodePtr _ps_node, ID2IDSet &_pairs ) {
 
231
void PS_find_probes_for_pairs(const PS_NodePtr& _ps_node, ID2IDSet &_pairs) {
224
232
    SpeciesID id         = _ps_node->getNum();
225
233
    bool      has_probes = _ps_node->hasProbes();
226
234
 
227
235
    //
228
236
    // append ID to path
229
237
    //
230
 
    __PATHSET->insert( id );
 
238
    __PATHSET->insert(id);
231
239
 
232
240
    //
233
 
    // dont look at path until ID is greater than lowest ID in the set of ID-pairs
 
241
    // don't look at path until ID is greater than lowest ID in the set of ID-pairs
234
242
    //
235
243
    if ((id >= __ONEMATCH_MIN_ID) && has_probes) {
236
244
        for (ID2IDSetCIter pair=_pairs.begin(); pair != _pairs.end(); ++pair) {
237
245
            // look for pair-IDs in the path
238
 
            bool found_first  = __PATHSET->find( pair->first  ) != __PATHSET->end();
239
 
            bool found_second = __PATHSET->find( pair->second ) != __PATHSET->end();
 
246
            bool found_first  = __PATHSET->find(pair->first) != __PATHSET->end();
 
247
            bool found_second = __PATHSET->find(pair->second) != __PATHSET->end();
240
248
            if (found_first ^ found_second) { // ^ is XOR
241
 
                printf( "found path for (%i,%i) at %p ", pair->first, pair->second,&(*_ps_node) );
 
249
                printf("found path for (%i,%i) at %p ", pair->first, pair->second, &(*_ps_node));
242
250
                _ps_node->printOnlyMe();
243
 
                (*__PAIR2PATH)[ *pair ] = *__PATHSET;   // store path
244
 
                _pairs.erase( pair );                   // remove found pair
245
 
                // scan pairs for new min,max IDs
246
 
                if ((pair->first == __ONEMATCH_MIN_ID) || (pair->second == __ONEMATCH_MAX_ID)) {
 
251
                (*__PAIR2PATH)[*pair] = *__PATHSET;     // store path
 
252
                
 
253
                bool scanMinMax = (pair->first == __ONEMATCH_MIN_ID) || (pair->second == __ONEMATCH_MAX_ID);
 
254
                _pairs.erase(pair);                     // remove found pair
 
255
                if (scanMinMax) {
 
256
                    // scan pairs for new min,max IDs
247
257
                    __ONEMATCH_MIN_ID = __MAX_ID;
248
258
                    __ONEMATCH_MAX_ID = -1;
249
259
                    for (ID2IDSetCIter p=_pairs.begin(); p != _pairs.end(); ++p) {
250
260
                        if (p->first  < __ONEMATCH_MIN_ID) __ONEMATCH_MIN_ID = p->first;
251
261
                        if (p->second > __ONEMATCH_MAX_ID) __ONEMATCH_MAX_ID = p->second;
252
262
                    }
253
 
                    printf( " new MIN,MAX (%d,%d)", __ONEMATCH_MIN_ID, __ONEMATCH_MAX_ID );
 
263
                    printf(" new MIN,MAX (%d,%d)", __ONEMATCH_MIN_ID, __ONEMATCH_MAX_ID);
254
264
                }
255
 
                printf( "\n" );
 
265
                printf("\n");
256
266
            }
257
267
        }
258
268
    }
259
 
        
 
269
 
260
270
    //
261
271
    // step down the children unless all paths are found
262
272
    // if either ID is lower than highest ID in the set of ID-pairs
266
276
        for (PS_NodeMapConstIterator i = _ps_node->getChildrenBegin();
267
277
             (i != _ps_node->getChildrenEnd()) && (!_pairs.empty());
268
278
             ++i) {
269
 
            PS_find_probes_for_pairs( i->second, _pairs );
 
279
            PS_find_probes_for_pairs(i->second, _pairs);
270
280
        }
271
281
    }
272
282
 
273
283
    //
274
284
    // remove ID from path
275
285
    //
276
 
    __PATHSET->erase( id );
 
286
    __PATHSET->erase(id);
277
287
}
278
288
 
279
 
void PS_print_and_evaluate_map( const PS_NodePtr _root_node, const char *_result_filename ) {
 
289
static void PS_print_and_evaluate_map(const PS_NodePtr& _root_node, const char *_result_filename) {
280
290
    //
281
291
    // print and evaluate bitmap
282
292
    //
283
 
    printf( "\n\n----------------- bitmap ---------------\n\n" );
 
293
    printf("\n\n----------------- bitmap ---------------\n\n");
284
294
    SpeciesID smaller_id;
285
295
    SpeciesID bigger_id;
286
296
    ID2IDSet  noMatch;
290
300
    __ONEMATCH_MIN_ID = __MAX_ID;
291
301
    __ONEMATCH_MAX_ID = __MIN_ID;
292
302
    for (SpeciesID id1 = __MIN_ID; id1 <= __MAX_ID; ++id1) {
293
 
//         printf( "[%6i] ",id1 );
294
303
        for (SpeciesID id2 = __MIN_ID; id2 <= id1; ++id2) {
295
304
            smaller_id = (id1 < id2) ? id1 : id2;
296
305
            bigger_id  = (id1 < id2) ? id2 : id1;
297
 
            bit1       = __MAP->get( smaller_id, bigger_id );
298
 
            bit2       = __MAP->get( bigger_id, smaller_id );
 
306
            bit1       = __MAP->get(smaller_id, bigger_id);
 
307
            bit2       = __MAP->get(bigger_id, smaller_id);
299
308
            if (bit1 && bit2) {
300
 
//                 printf( "2" );
301
 
            } else if (bit1) {
302
 
//                 printf( "1" );
303
 
                oneMatch.insert( ID2IDPair(smaller_id,bigger_id) );
 
309
            }
 
310
            else if (bit1) {
 
311
                oneMatch.insert(ID2IDPair(smaller_id, bigger_id));
304
312
                if (smaller_id < __ONEMATCH_MIN_ID) __ONEMATCH_MIN_ID = smaller_id;
305
313
                if (bigger_id  > __ONEMATCH_MAX_ID) __ONEMATCH_MAX_ID = bigger_id;
306
 
            } else {
307
 
//                 printf( "0" );
308
 
                if (id1 != id2) noMatch.insert( ID2IDPair(smaller_id,bigger_id) ); // there are no probes to distinguish a species from itself .. obviously
 
314
            }
 
315
            else {
 
316
                if (id1 != id2) noMatch.insert(ID2IDPair(smaller_id, bigger_id));  // there are no probes to distinguish a species from itself .. obviously
309
317
            }
310
318
        }
311
 
//         printf( "\n" );
312
319
    }
313
 
    printf( "(enter to continue)\n" );
314
 
//    getchar();
 
320
    printf("(enter to continue)\n");
315
321
 
316
 
    printf( "\n\n----------------- no matches ---------------\n\n" );
 
322
    printf("\n\n----------------- no matches ---------------\n\n");
317
323
    if (!_result_filename) {
318
324
        for (ID2IDSetCIter i = noMatch.begin(); i != noMatch.end(); ++i) {
319
 
            printf( "%6i %6i\n", i->first, i->second );
 
325
            printf("%6i %6i\n", i->first, i->second);
320
326
        }
321
327
    }
322
 
    printf( "%zu no matches\n(enter to continue)\n", noMatch.size() );
323
 
//    getchar();
 
328
    printf("%zu no matches\n(enter to continue)\n", noMatch.size());
324
329
 
325
 
    printf( "\n\n----------------- one match ---------------\n\n" );
 
330
    printf("\n\n----------------- one match ---------------\n\n");
326
331
    if (!_result_filename) {
327
332
        for (ID2IDSetCIter i = oneMatch.begin(); i != oneMatch.end(); ++i) {
328
 
            printf( "%6i %6i\n", i->first, i->second );
 
333
            printf("%6i %6i\n", i->first, i->second);
329
334
        }
330
335
    }
331
 
    printf( "%zu one matches\n(enter to continue)\n", oneMatch.size() );
332
 
//    getchar();
 
336
    printf("%zu one matches\n(enter to continue)\n", oneMatch.size());
333
337
    //
334
338
    // find paths for pairs
335
339
    //
338
342
    int c = 0;
339
343
    for (PS_NodeMapConstIterator i = _root_node->getChildrenBegin();
340
344
         (i != _root_node->getChildrenEnd()) && (!oneMatch.empty());
341
 
         ++i,++c ) {
342
 
        if ((c < 50) || (c % 100 == 0)) printf( "PS_find_probes_for_pairs( %i ) : %i of %zu\n", i->first, c+1, _root_node->countChildren() );
343
 
        PS_find_probes_for_pairs( i->second, oneMatch );
 
345
         ++i, ++c) {
 
346
        if ((c < 50) || (c % 100 == 0)) printf("PS_find_probes_for_pairs( %i ) : %i of %zu\n", i->first, c+1, _root_node->countChildren());
 
347
        PS_find_probes_for_pairs(i->second, oneMatch);
344
348
    }
345
349
    //
346
350
    // print paths
348
352
    for (IDID2IDSetMapCIter i = __PAIR2PATH->begin();
349
353
         i != __PAIR2PATH->end();
350
354
         ++i) {
351
 
        printf( "\nPair (%i,%i) Setsize (%zu)", i->first.first, i->first.second, i->second.size() );
 
355
        printf("\nPair (%i,%i) Setsize (%zu)", i->first.first, i->first.second, i->second.size());
352
356
        PS_NodePtr current_node = _root_node;
353
357
        long c2 = 0;
354
358
        for (IDSetCIter path_id=i->second.begin();
355
 
             path_id !=i->second.end();
356
 
             ++path_id,++c2) {
357
 
            current_node = current_node->getChild( *path_id ).second;
358
 
            if (c2 % 10 == 0) printf( "\n" );
359
 
            printf( "%6i%s ", *path_id, (current_node->hasProbes()) ? ((current_node->hasInverseProbes()) ? "*" : "+") : "-" );
 
359
             path_id != i->second.end();
 
360
             ++path_id, ++c2) {
 
361
            current_node = current_node->getChild(*path_id).second;
 
362
            if (c2 % 10 == 0) printf("\n");
 
363
            printf("%6i%s ", *path_id, (current_node->hasProbes()) ? ((current_node->hasInverseProbes()) ? "*" : "+") : "-");
360
364
        }
361
 
        printf( "\nFinal Node : %p ", &(*current_node) );
 
365
        printf("\nFinal Node : %p ", &(*current_node));
362
366
        current_node->printOnlyMe();
363
 
        printf( "\n" );
 
367
        printf("\n");
364
368
    }
365
 
    printf( "\n%zu paths\n", __PAIR2PATH->size() );
 
369
    printf("\n%zu paths\n", __PAIR2PATH->size());
366
370
    //
367
371
    // oups
368
372
    //
369
373
    if (!oneMatch.empty()) {
370
 
        printf( "did not find a path for these :\n" );
 
374
        printf("did not find a path for these :\n");
371
375
        for (ID2IDSetCIter i = oneMatch.begin(); i != oneMatch.end(); ++i) {
372
 
            printf( "%6i %6i\n", i->first, i->second );
 
376
            printf("%6i %6i\n", i->first, i->second);
373
377
        }
374
378
    }
375
379
    //
376
380
    // make preset map
377
381
    //
378
 
    PS_BitMap_Counted *preset = new PS_BitMap_Counted( false, __MAX_ID+1 );
 
382
    PS_BitMap_Counted *preset = new PS_BitMap_Counted(false, __MAX_ID+1);
379
383
    // set bits for no matches
380
384
    for (ID2IDSetCIter pair=noMatch.begin(); pair != noMatch.end(); ++pair) {
381
 
        preset->setTrue( pair->second, pair->first );
 
385
        preset->setTrue(pair->second, pair->first);
382
386
    }
383
387
    // iterate over paths
384
388
    for (IDID2IDSetMapCIter i = __PAIR2PATH->begin();
397
401
            for (IDSetCIter path_id = i->second.begin(); path_id != i->second.end(); ++path_id) {
398
402
                if (id == *path_id) continue;   // obviously a probe cant differ a species from itself
399
403
                if (id > *path_id) {
400
 
                    preset->setTrue( id,*path_id );
401
 
                } else {
402
 
                    preset->setTrue( *path_id,id );
 
404
                    preset->setTrue(id, *path_id);
 
405
                }
 
406
                else {
 
407
                    preset->setTrue(*path_id, id);
403
408
                }
404
409
            }
405
410
        }
410
415
    // save results
411
416
    //
412
417
    if (_result_filename) {
413
 
        PS_FileBuffer *result_file = new PS_FileBuffer( _result_filename, PS_FileBuffer::WRITEONLY );
 
418
        PS_FileBuffer *result_file = new PS_FileBuffer(_result_filename, PS_FileBuffer::WRITEONLY);
414
419
        // no matches
415
 
        printf( "saving no matches to %s...\n", _result_filename );
416
 
        result_file->put_long( noMatch.size() );
 
420
        printf("saving no matches to %s...\n", _result_filename);
 
421
        result_file->put_long(noMatch.size());
417
422
        for (ID2IDSetCIter i = noMatch.begin(); i != noMatch.end(); ++i) {
418
 
            result_file->put_int( i->first );
419
 
            result_file->put_int( i->second );
 
423
            result_file->put_int(i->first);
 
424
            result_file->put_int(i->second);
420
425
        }
421
426
        // one matches
422
 
        printf( "saving one matches to %s...\n", _result_filename );
423
 
        result_file->put_long( __PAIR2PATH->size() );
 
427
        printf("saving one matches to %s...\n", _result_filename);
 
428
        result_file->put_long(__PAIR2PATH->size());
424
429
        for (IDID2IDSetMapCIter i = __PAIR2PATH->begin(); i != __PAIR2PATH->end(); ++i) {
425
 
            result_file->put_int( i->first.first );
426
 
            result_file->put_int( i->first.second );
427
 
            result_file->put_long( i->second.size() );
428
 
            for (IDSetCIter path_id=i->second.begin(); path_id !=i->second.end(); ++path_id) {
429
 
                result_file->put_int( *path_id );
 
430
            result_file->put_int(i->first.first);
 
431
            result_file->put_int(i->first.second);
 
432
            result_file->put_long(i->second.size());
 
433
            for (IDSetCIter path_id=i->second.begin(); path_id != i->second.end(); ++path_id) {
 
434
                result_file->put_int(*path_id);
430
435
            }
431
436
        }
432
437
        // preset bitmap
433
 
        printf( "saving preset bitmap to %s...\n", _result_filename );
434
 
        preset->save( result_file );
 
438
        printf("saving preset bitmap to %s...\n", _result_filename);
 
439
        preset->save(result_file);
435
440
        delete result_file;
436
441
    }
437
442
    delete preset;
442
447
//  ====================================================
443
448
//  ====================================================
444
449
 
445
 
int main( int argc,  char *argv[] ) {
 
450
int main(int argc,   char *argv[]) {
446
451
 
447
452
   // open probe-set-database
448
453
    if (argc < 2) {
449
 
        printf( "Missing argument\n Usage %s <database name> [[-]bitmap filename] [+result filename]\n ", argv[0] );
450
 
        printf( "if bitmap filename begins with '-' it is loaded, else its created\n " );
451
 
        printf( "result filename MUST be preceded by '+'\n" );
452
 
        exit( 1 );
 
454
        printf("Missing argument\n Usage %s <database name> [[-]bitmap filename] [+result filename]\n ", argv[0]);
 
455
        printf("if bitmap filename begins with '-' it is loaded, else its created\n ");
 
456
        printf("result filename MUST be preceded by '+'\n");
 
457
        exit(1);
453
458
    }
454
459
 
455
460
    const char *input_DB_name   = argv[1];
459
464
    if (argc > 2) {
460
465
        if (argv[2][0] == '+') {
461
466
            result_filename = argv[2]+1;
462
 
        } else {
 
467
        }
 
468
        else {
463
469
            bitmap_filename = argv[2];
464
470
        }
465
471
    }
466
472
    if (argc > 3) {
467
473
        if (argv[3][0] == '+') {
468
474
            result_filename = argv[3]+1;
469
 
        } else {
 
475
        }
 
476
        else {
470
477
            bitmap_filename = argv[3];
471
478
        }
472
479
    }
473
 
    
 
480
 
474
481
    struct tms before;
475
 
    times( &before );
476
 
    printf( "Opening probe-set-database '%s'..\n", input_DB_name );
477
 
    PS_Database *db = new PS_Database( input_DB_name, PS_Database::READONLY );
 
482
    times(&before);
 
483
    printf("Opening probe-set-database '%s'..\n", input_DB_name);
 
484
    PS_Database *db = new PS_Database(input_DB_name, PS_Database::READONLY);
478
485
    db->load();
479
486
    __MAX_ID = db->getMaxID();
480
487
    __MIN_ID = db->getMinID();
481
 
    PS_print_time_diff( &before, "(enter to continue)  " );
482
 
//    getchar();
 
488
    PS_print_time_diff(&before, "(enter to continue)  ");
483
489
 
484
 
    __MAP = new PS_BitMap_Fast( false, __MAX_ID+1 );
 
490
    __MAP = new PS_BitMap_Fast(false, __MAX_ID+1);
485
491
    if (!bitmap_filename || (bitmap_filename[0] != '-')) {
486
 
        printf( "detecting..\n" ); fflush( stdout );
487
 
        __NODES_LEFT = (char*)malloc( 61 );
488
 
        memset( __NODES_LEFT, ' ', 60 );
489
 
        __NODES_LEFT[ 60 ] = '\x00';
490
 
        PS_detect_weak_differences( db->getRootNode() );
491
 
        free( __NODES_LEFT );
 
492
        printf("detecting..\n"); fflush(stdout);
 
493
        __NODES_LEFT = (char*)malloc(61);
 
494
        memset(__NODES_LEFT, ' ', 60);
 
495
        __NODES_LEFT[60] = '\x00';
 
496
        PS_detect_weak_differences(db->getRootNode());
 
497
        free(__NODES_LEFT);
492
498
        if (bitmap_filename) {
493
 
            printf( "saving bitmap to file %s\n", bitmap_filename );
494
 
            PS_FileBuffer *mapfile = new PS_FileBuffer( bitmap_filename, PS_FileBuffer::WRITEONLY );
495
 
            __MAP->save( mapfile );
 
499
            printf("saving bitmap to file %s\n", bitmap_filename);
 
500
            PS_FileBuffer *mapfile = new PS_FileBuffer(bitmap_filename, PS_FileBuffer::WRITEONLY);
 
501
            __MAP->save(mapfile);
496
502
            delete mapfile;
497
503
        }
498
 
    } else if (bitmap_filename) {
499
 
        printf( "loading bitmap from file %s\n",bitmap_filename+1 );
500
 
        PS_FileBuffer *mapfile = new PS_FileBuffer( bitmap_filename+1, PS_FileBuffer::READONLY );
501
 
        __MAP->load( mapfile );
 
504
    }
 
505
    else if (bitmap_filename) {
 
506
        printf("loading bitmap from file %s\n", bitmap_filename+1);
 
507
        PS_FileBuffer *mapfile = new PS_FileBuffer(bitmap_filename+1, PS_FileBuffer::READONLY);
 
508
        __MAP->load(mapfile);
502
509
        delete mapfile;
503
510
    }
504
 
    printf( "(enter to continue)\n" );
505
 
//    getchar();
 
511
    printf("(enter to continue)\n");
506
512
 
507
 
    times( &before );
508
 
    PS_print_and_evaluate_map( db->getRootNode(), result_filename );
509
 
    PS_print_time_diff( &before, "(enter to continue)  " );
510
 
//    getchar();
 
513
    times(&before);
 
514
    PS_print_and_evaluate_map(db->getRootNode(), result_filename);
 
515
    PS_print_time_diff(&before, "(enter to continue)  ");
511
516
    delete __MAP;
512
 
    
513
 
    printf( "removing database from memory\n" );
 
517
 
 
518
    printf("removing database from memory\n");
514
519
    delete db;
515
 
    printf( "(enter to continue)\n" );
516
 
//    getchar();
 
520
    printf("(enter to continue)\n");
517
521
 
518
522
    return 0;
519
523
}