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