1
2 /******************************************************************************
3 *
4 * This file is part of meryl-utility, a collection of miscellaneous code
5 * used by Meryl, Canu and others.
6 *
7 * This software is based on:
8 * 'Canu' v2.0 (https://github.com/marbl/canu)
9 * which is based on:
10 * 'Celera Assembler' r4587 (http://wgs-assembler.sourceforge.net)
11 * the 'kmer package' r1994 (http://kmer.sourceforge.net)
12 *
13 * Except as indicated otherwise, this is a 'United States Government Work',
14 * and is released in the public domain.
15 *
16 * File 'README.licenses' in the root directory of this distribution
17 * contains full conditions and disclaimers.
18 */
19
20 #include "runtime.H"
21 #include "intervalList.H"
22
23 #include "mt19937ar.H"
24
25 // g++ -o intervalListTest -I.. -I. intervalListTest.C
26
27 int
main(int argc,char ** argv)28 main(int argc, char **argv) {
29
30 AS_configure(argc, argv);
31
32 // Boring basic test.
33 if (0) {
34 intervalList<int32> t1;
35
36 t1.add(0, 10);
37 t1.add(11,7);
38 t1.add(20, 8);
39
40 fprintf(stderr, "BEFORE:\n");
41 for (uint32 ii=0; ii<t1.numberOfIntervals(); ii++)
42 fprintf(stderr, "%2d %3d-%3d\n", ii, t1.lo(ii), t1.hi(ii));
43
44 t1.merge(-1);
45
46 fprintf(stderr, "AFTER:\n");
47 for (uint32 ii=0; ii<t1.numberOfIntervals(); ii++)
48 fprintf(stderr, "%2d %3d-%3d\n", ii, t1.lo(ii), t1.hi(ii));
49 }
50
51
52 if (1) {
53 mtRandom mt(strtouint32(argv[1]));
54
55 // About 6.5 minutes per million, so this should be about an hour.
56 uint32 iterMax = 935000;
57
58 for (uint32 iter=0; iter<iterMax; iter++) {
59 uint32 numIntervals = mt.mtRandom32() % 5000;
60 uint32 maxLen = 1 + mt.mtRandom32() % 1000;
61 uint32 maxBgn = 1 + mt.mtRandom32() % 50000;
62 uint32 *depth = new uint32 [maxBgn + maxLen];
63
64 memset(depth, 0, sizeof(uint32) * (maxBgn + maxLen));
65
66 if (iter % 1000 == 0)
67 fprintf(stderr, "%10u/%10u: %3u intervals, each up to %4u long, coords up to %4u\n",
68 iter, iterMax,
69 numIntervals, maxLen, maxBgn);
70
71 intervalList<uint32> il;
72
73 // Add intervals to the list.
74 // Sum depths explicitly.
75 for (uint32 ii=0; ii<numIntervals; ii++) {
76 uint32 bgn = mt.mtRandom32() % maxBgn; // bgn between 0 and maxBgn
77 uint32 len = mt.mtRandom32() % maxLen; // len between 0 and maxLen
78
79 il.add(bgn, len);
80
81 for (uint32 xx=bgn; xx<bgn+len; xx++)
82 depth[xx]++;
83
84 //fprintf(stderr, "IL %u - %u\n", bgn, bgn+len);
85 }
86
87 // Convert intervals to depths.
88 intervalDepth<uint32> de(il);
89
90 // Over all the depth regions, subtract the computed depth from
91 // the explicit depth.
92 for (uint32 xx=0; xx<de.numberOfIntervals(); xx++) {
93 uint32 bgn = de.lo(xx);
94 uint32 end = de.hi(xx);
95 uint32 dpt = de.depth(xx);
96
97 //fprintf(stderr, "ID %u - %u depth %u\n", bgn, end, dpt);
98
99 for (uint32 cc=bgn; cc<end; cc++) {
100 //if (cc < 30)
101 // fprintf(stderr, "depth[%u] = %u -> %u\n", cc, depth[cc], depth[cc] - dpt);
102 depth[cc] -= dpt;
103 }
104 }
105
106 // Every explicit depth should now be zero, even the ones
107 // not covered.
108 for (uint32 cc=0; cc<maxBgn + maxLen; cc++) {
109 if (depth[cc] != 0)
110 fprintf(stderr, "ERROR: depth[%u] = %u in iter %u\n", cc, depth[cc], iter);
111 assert(depth[cc] == 0);
112 }
113
114 delete [] depth;
115 }
116 }
117
118 exit(0);
119 }
120