1 #include <gtest/gtest.h>
2 #include <unit_tests/test_main.hpp>
3 #include <stdexcept>
4 #include <stdlib.h>
5 #include <jellyfish/rectangular_binary_matrix.hpp>
6 #include <jellyfish/mer_dna.hpp>
7
8 #ifdef HAVE_CONFIG_H
9 #include <config.h>
10 #endif
11
12 #ifdef HAVE_INT128
13 #include <jellyfish/int128.hpp>
14 #endif
15
16 namespace {
17 using jellyfish::RectangularBinaryMatrix;
18 using jellyfish::mer_dna;
19
allocate_matrix(unsigned int r,unsigned c)20 static bool allocate_matrix(unsigned int r, unsigned c) {
21 RectangularBinaryMatrix m(r, c);
22 return m.is_zero();
23 }
24
TEST(RectangularBinaryMatrix,InitSizes)25 TEST(RectangularBinaryMatrix, InitSizes) {
26 RectangularBinaryMatrix m(5, 60);
27
28 EXPECT_EQ(5u, m.r());
29 EXPECT_EQ(60u, m.c());
30 EXPECT_TRUE(m.is_zero());
31
32 EXPECT_THROW(allocate_matrix(100, 100), std::out_of_range);
33 EXPECT_THROW(allocate_matrix(0, 100), std::out_of_range);
34 EXPECT_THROW(allocate_matrix(10, 0), std::out_of_range);
35 }
36
TEST(RectangularBinaryMatrix,Copy)37 TEST(RectangularBinaryMatrix, Copy) {
38 RectangularBinaryMatrix m1(5, 55);
39 m1.randomize(random_bits);
40
41 RectangularBinaryMatrix m2(m1);
42 RectangularBinaryMatrix m3(6, 66);
43 RectangularBinaryMatrix m4(5, 55);
44
45 EXPECT_TRUE(!m1.is_zero());
46 EXPECT_TRUE(m1 == m2);
47 EXPECT_TRUE(!(m1 == m3));
48 EXPECT_TRUE(!(m1 == m4));
49 m4 = m1;
50 EXPECT_TRUE(m1 == m4);
51 }
52
TEST(RectangularBinaryMatrix,InitRaw)53 TEST(RectangularBinaryMatrix, InitRaw) {
54 const unsigned int nb_col = 80;
55 uint64_t raw[nb_col];
56 for(unsigned int i = 0; i < nb_col; ++i)
57 raw[i] = random_bits();
58 const RectangularBinaryMatrix m(raw, 19, nb_col);
59 EXPECT_EQ(19u, m.r());
60 EXPECT_EQ(80u, m.c());
61 const uint64_t mask = ((uint64_t)1 << 19) - 1;
62 for(unsigned int i = 0; i < nb_col; ++i)
63 EXPECT_EQ(raw[i] & mask, m[i]);
64 }
65
TEST(RectangularBinaryMatrix,LowIdentity)66 TEST(RectangularBinaryMatrix, LowIdentity) {
67 for(int r = 2; r < 64; r += 2) {
68 for(int c = 2; c < 100; c += 2) {
69 SCOPED_TRACE(::testing::Message() << "matrix " << r << "x" << c);
70 RectangularBinaryMatrix m(r, c); // Matrix should be zeroed out
71 mer_dna::k(c);
72 mer_dna v;
73 v.randomize();
74
75 EXPECT_FALSE(m.is_low_identity());
76 m.init_low_identity();
77 EXPECT_TRUE(m.is_low_identity());
78
79 uint64_t res = m.times(v);
80 EXPECT_EQ(v.get_bits(0, std::min(r, c)), res);
81
82 RectangularBinaryMatrix m2 = RectangularBinaryMatrix::identity(r);
83 uint64_t res2 = m2.times(v);
84 EXPECT_EQ(v.get_bits(0, r), res2);
85 }
86 }
87 }
88
89 /******************************
90 * Matrix Vector product
91 ******************************/
92 class MatrixVectorProd : public ::testing::TestWithParam< ::std::tr1::tuple<int, int> > {
93 public:
94 unsigned int row, col;
95 RectangularBinaryMatrix m;
MatrixVectorProd()96 MatrixVectorProd() :
97 row(::std::tr1::get<0>(GetParam())),
98 col(::std::tr1::get<1>(GetParam())),
99 m(row, col)
100 {
101 m.randomize(random_bits);
102 }
103 };
104
TEST_P(MatrixVectorProd,Checks)105 TEST_P(MatrixVectorProd, Checks) {
106 EXPECT_EQ(row, m.r());
107 EXPECT_EQ(col, m.c());
108 }
109
TEST_P(MatrixVectorProd,AllOnes)110 TEST_P(MatrixVectorProd, AllOnes) {
111 uint64_t v[2], res = 0;
112 v[0] = v[1] = (uint64_t)-1;
113 for(unsigned int i = 0; i < m.c(); ++i)
114 res ^= m[i];
115 EXPECT_EQ(res, m.times_loop(v));
116 #ifdef HAVE_INT128
117 EXPECT_EQ(res, m.times_128(v));
118 #endif
119 #ifdef HAVE_SSE
120 EXPECT_EQ(res, m.times_sse(v));
121 #endif
122 }
123
TEST_P(MatrixVectorProd,EveryOtherOnes)124 TEST_P(MatrixVectorProd, EveryOtherOnes) {
125 uint64_t v[2], res = 0;
126 v[0] = 0xaaaaaaaaaaaaaaaaUL;
127 v[1] = 0xaaaaaaaaaaaaaaaaUL;
128 for(unsigned int i = 0; i < m.c(); i += 2)
129 res ^= m[i];
130 EXPECT_EQ(res, m.times_loop(v));
131 #ifdef HAVE_INT128
132 EXPECT_EQ(res, m.times_128(v));
133 #endif
134 #ifdef HAVE_SSE
135 EXPECT_EQ(res, m.times_sse(v));
136 #endif
137
138 v[0] >>= 1;
139 v[1] >>= 1;
140 res = 0;
141 for(unsigned int i = 1; i < m.c(); i += 2)
142 res ^= m[i];
143 EXPECT_EQ(res, m.times_loop(v));
144 #ifdef HAVE_INT128
145 EXPECT_EQ(res, m.times_128(v));
146 #endif
147 #ifdef HAVE_SSE
148 EXPECT_EQ(res, m.times_sse(v));
149 #endif
150 }
151
152 #if HAVE_SSE || HAVE_INT128
TEST_P(MatrixVectorProd,Optimizations)153 TEST_P(MatrixVectorProd, Optimizations) {
154 static const int nb_tests = 100;
155 const unsigned int nb_words = col / 64 + (col % 64 != 0);
156 uint64_t v[nb_words];
157
158 for(int i = 0; i < nb_tests; ++i) {
159 // unsigned int r = 2 * (random() % 31 + 1);
160 // unsigned int c = 2 * (random() % 100) + r;
161
162 // RectangularBinaryMatrix m(r, c);
163 // m.randomize(random_bits);
164
165 for(unsigned int j = 0; j < nb_words; ++j)
166 v[j] = random_bits();
167
168 uint64_t res = m.times_loop((uint64_t*)v);
169 #ifdef HAVE_SSE
170 EXPECT_EQ(res, m.times_sse((uint64_t*)v));
171 #endif
172 #ifdef HAVE_INT128
173 EXPECT_EQ(res, m.times_128((uint64_t*)v));
174 #endif
175 }
176 }
177 #endif // HAVE_SSE || HAVE_INT128
178
179 INSTANTIATE_TEST_CASE_P(MatrixVectorProdTest, MatrixVectorProd, ::testing::Combine(::testing::Range(1, 65, 1), // rows
180 ::testing::Range(2, 100, 2))); // cols
181
182 /******************************
183 * Matrix product and inverse
184 ******************************/
TEST(PseudoProduct,Dimensions)185 TEST(PseudoProduct, Dimensions) {
186 RectangularBinaryMatrix m(30, 100), m1(32, 100), m2(30, 98);
187
188 EXPECT_THROW(m.pseudo_multiplication(m1), std::domain_error);
189 EXPECT_THROW(m.pseudo_multiplication(m2), std::domain_error);
190 }
191
TEST(PseudoProduct,Identity)192 TEST(PseudoProduct, Identity) {
193 RectangularBinaryMatrix m(30, 100), i(30, 100);
194 i.init_low_identity();
195 m.randomize(random_bits);
196
197 EXPECT_TRUE(i.pseudo_multiplication(m) == m);
198 }
199
TEST(PseudoProduct,Parity)200 TEST(PseudoProduct, Parity) {
201 const unsigned int col_sizes[6] = { 50, 70, 126, 130, 64, 128 };
202 const unsigned int nb_rows = 30;
203
204 for(unsigned int k = 0; k < sizeof(col_sizes) / sizeof(unsigned int); ++k) {
205 const unsigned int nb_cols = col_sizes[k];
206 uint64_t *cols = new uint64_t[nb_cols];
207 RectangularBinaryMatrix p(nb_rows, nb_cols);
208
209 for(unsigned int j = 18; j < 19; ++j) {
210 const uint64_t bits = ((uint64_t)1 << j) - 1;
211 unsigned int i;
212 for(i = 0; i < nb_cols; ++i)
213 cols[i] = bits;
214 RectangularBinaryMatrix m(cols, nb_rows, nb_cols);
215
216 p = m.pseudo_multiplication(m);
217 for(i = 0; i < nb_cols - nb_rows; ++i)
218 EXPECT_EQ(__builtin_parity(bits) ? (uint64_t)0 : bits, p[i]);
219 for( ; i < nb_cols; ++i)
220 EXPECT_EQ(__builtin_parity(bits) ? bits : (uint64_t)0, p[i]);
221 }
222 delete [] cols;
223 }
224 }
225
TEST(PseudoProduct,Inverse)226 TEST(PseudoProduct, Inverse) {
227 int full_rank = 0, singular = 0;
228 for(unsigned int i = 0; i < 500; ++i) {
229 unsigned int r = random() % 63 + 1;
230 unsigned int c = 2 * ((random() % 100) + 1);
231 SCOPED_TRACE(::testing::Message() << "Dimension " << r << "x" << c);
232 RectangularBinaryMatrix m(r, c);
233 m.randomize(random_bits);
234 RectangularBinaryMatrix s(m);
235 unsigned int rank = m.pseudo_rank();
236 if(rank != c) {
237 ++singular;
238 EXPECT_THROW(m.pseudo_inverse(), std::domain_error);
239 } else {
240 ++full_rank;
241 RectangularBinaryMatrix inv(m);
242 EXPECT_NO_THROW(inv = m.pseudo_inverse());
243 RectangularBinaryMatrix i = inv.pseudo_multiplication(m);
244 EXPECT_TRUE(i.is_low_identity());
245 }
246 EXPECT_TRUE(s == m);
247 }
248 EXPECT_EQ(500, full_rank + singular);
249 EXPECT_NE(0, full_rank);
250 }
251
TEST(PseudoProduct,Rank)252 TEST(PseudoProduct, Rank) {
253 RectangularBinaryMatrix m(50, 100);
254 for(unsigned int i = 0; i < 10; ++i) {
255 m.randomize(random_bits);
256 RectangularBinaryMatrix s(m);
257 unsigned int rank = m.pseudo_rank();
258 EXPECT_TRUE(rank <= m.c());
259 EXPECT_TRUE(s == m);
260 }
261 }
262
TEST(PseudoProduct,InitRandom)263 TEST(PseudoProduct, InitRandom) {
264 RectangularBinaryMatrix m(50, 100);
265 for(unsigned int i = 0; i < 10; ++i) {
266 RectangularBinaryMatrix im(m.randomize_pseudo_inverse(random_bits));
267 EXPECT_EQ((unsigned int)m.c(), m.pseudo_rank());
268 EXPECT_EQ((unsigned int)m.c(), im.pseudo_rank());
269 EXPECT_TRUE((m.pseudo_multiplication(im)).is_low_identity());
270 }
271 }
272
TEST(PseudoProduct,VectorInverseMultiplication)273 TEST(PseudoProduct, VectorInverseMultiplication) {
274 RectangularBinaryMatrix m(50, 132);
275 RectangularBinaryMatrix im(m.randomize_pseudo_inverse());
276 RectangularBinaryMatrix cim(im);
277 EXPECT_TRUE(cim.pseudo_multiplication(m).is_low_identity());
278
279 mer_dna::k(66);
280 mer_dna v, iv, iv2;
281 EXPECT_EQ(m.nb_words(), v.nb_words());
282
283 for(int i = 0; i < 100; ++i) {
284 SCOPED_TRACE(::testing::Message() << "i=" << i);
285 v.randomize();
286 uint64_t hash = m.times(v);
287 iv = v; iv.set_bits(0, 50, hash);
288
289 uint64_t lower = im.times(iv);
290 EXPECT_EQ(v.get_bits(0, 50), lower);
291 }
292 }
293
294
295 static const int speed_loop = 100000000;
TEST(MatrixProductSpeed,Loop)296 TEST(MatrixProductSpeed, Loop) {
297 RectangularBinaryMatrix m(50, 100);
298 const unsigned int nb_words = m.c() / 64 + (m.c() % 64 != 0);
299 uint64_t v[nb_words];
300 for(unsigned int j = 0; j < nb_words; ++j)
301 v[j] = random_bits();
302
303 volatile uint64_t res = 0;
304 for(int i = 0; i < speed_loop; ++i)
305 res ^= m.times_loop((uint64_t*)v);
306 }
307
308 #ifdef HAVE_SSE
TEST(MatrixProductSpeed,SSE)309 TEST(MatrixProductSpeed, SSE) {
310 RectangularBinaryMatrix m(50, 100);
311 const unsigned int nb_words = m.c() / 64 + (m.c() % 64 != 0);
312 uint64_t v[nb_words];
313 for(unsigned int j = 0; j < nb_words; ++j)
314 v[j] = random_bits();
315
316 volatile uint64_t res = 0;
317 for(int i = 0; i < speed_loop; ++i)
318 res ^= m.times_sse((uint64_t*)v);
319 }
320 #endif
321
322 #ifdef HAVE_INT128
TEST(MatrixProductSpeed,U128)323 TEST(MatrixProductSpeed, U128) {
324 RectangularBinaryMatrix m(50, 100);
325 const unsigned int nb_words = m.c() / 64 + (m.c() % 64 != 0);
326 uint64_t v[nb_words];
327 for(unsigned int j = 0; j < nb_words; ++j)
328 v[j] = random_bits();
329
330 volatile uint64_t res = 0;
331 for(int i = 0; i < speed_loop; ++i)
332 res ^= m.times_128((uint64_t*)v);
333 }
334 #endif
335
336 }
337