1 // Copyright (c) 2018-2019 Cem Bassoy
2 //
3 // Distributed under the Boost Software License, Version 1.0. (See
4 // accompanying file LICENSE_1_0.txt or copy at
5 // http://www.boost.org/LICENSE_1_0.txt)
6 //
7 // The authors gratefully acknowledge the support of
8 // Fraunhofer and Google in producing this work
9 // which started as a Google Summer of Code project.
10 //
11 // And we acknowledge the support from all contributors.
12
13
14 #include <iostream>
15 #include <algorithm>
16 #include <boost/numeric/ublas/tensor.hpp>
17 #include <boost/numeric/ublas/matrix.hpp>
18 #include <boost/numeric/ublas/vector.hpp>
19
20 #include <boost/test/unit_test.hpp>
21
22 #include "utility.hpp"
23
24 BOOST_AUTO_TEST_SUITE ( test_tensor_functions, * boost::unit_test::depends_on("test_tensor_contraction") )
25
26
27 using test_types = zip<int,long,float,double,std::complex<float>>::with_t<boost::numeric::ublas::first_order, boost::numeric::ublas::last_order>;
28
29 //using test_types = zip<int>::with_t<boost::numeric::ublas::first_order>;
30
31
32 struct fixture
33 {
34 using extents_type = boost::numeric::ublas::shape;
fixturefixture35 fixture()
36 : extents {
37 extents_type{1,1}, // 1
38 extents_type{1,2}, // 2
39 extents_type{2,1}, // 3
40 extents_type{2,3}, // 4
41 extents_type{2,3,1}, // 5
42 extents_type{4,1,3}, // 6
43 extents_type{1,2,3}, // 7
44 extents_type{4,2,3}, // 8
45 extents_type{4,2,3,5}} // 9
46 {
47 }
48 std::vector<extents_type> extents;
49 };
50
51
52
53
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_prod_vector,value,test_types,fixture)54 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_vector, value, test_types, fixture )
55 {
56 using namespace boost::numeric;
57 using value_type = typename value::first_type;
58 using layout_type = typename value::second_type;
59 using tensor_type = ublas::tensor<value_type,layout_type>;
60 using vector_type = typename tensor_type::vector_type;
61
62
63 for(auto const& n : extents){
64
65 auto a = tensor_type(n, value_type{2});
66
67 for(auto m = 0u; m < n.size(); ++m){
68
69 auto b = vector_type (n[m], value_type{1} );
70
71 auto c = ublas::prod(a, b, m+1);
72
73 for(auto i = 0u; i < c.size(); ++i)
74 BOOST_CHECK_EQUAL( c[i] , value_type(n[m]) * a[i] );
75
76 }
77 }
78 }
79
80
81
82
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_prod_matrix,value,test_types,fixture)83 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_matrix, value, test_types, fixture )
84 {
85 using namespace boost::numeric;
86 using value_type = typename value::first_type;
87 using layout_type = typename value::second_type;
88 using tensor_type = ublas::tensor<value_type,layout_type>;
89 using matrix_type = typename tensor_type::matrix_type;
90
91
92 for(auto const& n : extents) {
93
94 auto a = tensor_type(n, value_type{2});
95
96 for(auto m = 0u; m < n.size(); ++m){
97
98 auto b = matrix_type ( n[m], n[m], value_type{1} );
99
100 auto c = ublas::prod(a, b, m+1);
101
102 for(auto i = 0u; i < c.size(); ++i)
103 BOOST_CHECK_EQUAL( c[i] , value_type(n[m]) * a[i] );
104
105 }
106 }
107 }
108
109
110
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_prod_tensor_1,value,test_types,fixture)111 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_tensor_1, value, test_types, fixture )
112 {
113 using namespace boost::numeric;
114 using value_type = typename value::first_type;
115 using layout_type = typename value::second_type;
116 using tensor_type = ublas::tensor<value_type,layout_type>;
117
118 // left-hand and right-hand side have the
119 // the same number of elements
120
121 for(auto const& na : extents) {
122
123 auto a = tensor_type( na, value_type{2} );
124 auto b = tensor_type( na, value_type{3} );
125
126 auto const pa = a.rank();
127
128 // the number of contractions is changed.
129 for( auto q = 0ul; q <= pa; ++q) { // pa
130
131 auto phi = std::vector<std::size_t> ( q );
132
133 std::iota(phi.begin(), phi.end(), 1ul);
134
135 auto c = ublas::prod(a, b, phi);
136
137 auto acc = value_type(1);
138 for(auto i = 0ul; i < q; ++i)
139 acc *= a.extents().at(phi.at(i)-1);
140
141 for(auto i = 0ul; i < c.size(); ++i)
142 BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
143
144 }
145 }
146 }
147
148
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_prod_tensor_2,value,test_types,fixture)149 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_tensor_2, value, test_types, fixture )
150 {
151 using namespace boost::numeric;
152 using value_type = typename value::first_type;
153 using layout_type = typename value::second_type;
154 using tensor_type = ublas::tensor<value_type,layout_type>;
155
156
157 auto compute_factorial = [](auto const& p){
158 auto f = 1ul;
159 for(auto i = 1u; i <= p; ++i)
160 f *= i;
161 return f;
162 };
163
164 auto permute_extents = [](auto const& pi, auto const& na){
165 auto nb = na;
166 assert(pi.size() == na.size());
167 for(auto j = 0u; j < pi.size(); ++j)
168 nb[pi[j]-1] = na[j];
169 return nb;
170 };
171
172
173 // left-hand and right-hand side have the
174 // the same number of elements
175
176 for(auto const& na : extents) {
177
178 auto a = tensor_type( na, value_type{2} );
179 auto const pa = a.rank();
180
181
182 auto pi = std::vector<std::size_t>(pa);
183 auto fac = compute_factorial(pa);
184 std::iota( pi.begin(), pi.end(), 1 );
185
186 for(auto f = 0ul; f < fac; ++f)
187 {
188 auto nb = permute_extents( pi, na );
189 auto b = tensor_type( nb, value_type{3} );
190
191 // the number of contractions is changed.
192 for( auto q = 0ul; q <= pa; ++q) { // pa
193
194 auto phia = std::vector<std::size_t> ( q ); // concatenation for a
195 auto phib = std::vector<std::size_t> ( q ); // concatenation for b
196
197 std::iota(phia.begin(), phia.end(), 1ul);
198 std::transform( phia.begin(), phia.end(), phib.begin(),
199 [&pi] ( std::size_t i ) { return pi.at(i-1); } );
200
201 auto c = ublas::prod(a, b, phia, phib);
202
203 auto acc = value_type(1);
204 for(auto i = 0ul; i < q; ++i)
205 acc *= a.extents().at(phia.at(i)-1);
206
207 for(auto i = 0ul; i < c.size(); ++i)
208 BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
209
210 }
211
212 std::next_permutation(pi.begin(), pi.end());
213 }
214 }
215 }
216
217
218
219
220
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_inner_prod,value,test_types,fixture)221 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_inner_prod, value, test_types, fixture )
222 {
223 using namespace boost::numeric;
224 using value_type = typename value::first_type;
225 using layout_type = typename value::second_type;
226 using tensor_type = ublas::tensor<value_type,layout_type>;
227
228
229 for(auto const& n : extents) {
230
231 auto a = tensor_type(n, value_type(2));
232 auto b = tensor_type(n, value_type(1));
233
234 auto c = ublas::inner_prod(a, b);
235 auto r = std::inner_product(a.begin(),a.end(), b.begin(),value_type(0));
236
237 BOOST_CHECK_EQUAL( c , r );
238
239 }
240 }
241
242
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_norm,value,test_types,fixture)243 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_norm, value, test_types, fixture )
244 {
245 using namespace boost::numeric;
246 using value_type = typename value::first_type;
247 using layout_type = typename value::second_type;
248 using tensor_type = ublas::tensor<value_type,layout_type>;
249
250
251 for(auto const& n : extents) {
252
253 auto a = tensor_type(n);
254
255 auto one = value_type(1);
256 auto v = one;
257 for(auto& aa: a)
258 aa = v, v += one;
259
260
261 auto c = ublas::inner_prod(a, a);
262 auto r = std::inner_product(a.begin(),a.end(), a.begin(),value_type(0));
263
264 auto r2 = ublas::norm( (a+a) / 2 );
265
266 BOOST_CHECK_EQUAL( c , r );
267 BOOST_CHECK_EQUAL( std::sqrt( c ) , r2 );
268
269 }
270 }
271
272
BOOST_FIXTURE_TEST_CASE(test_tensor_real_imag_conj,fixture)273 BOOST_FIXTURE_TEST_CASE( test_tensor_real_imag_conj, fixture )
274 {
275 using namespace boost::numeric;
276 using value_type = float;
277 using complex_type = std::complex<value_type>;
278 using layout_type = ublas::first_order;
279
280 using tensor_complex_type = ublas::tensor<complex_type,layout_type>;
281 using tensor_type = ublas::tensor<value_type,layout_type>;
282
283 for(auto const& n : extents) {
284
285 auto a = tensor_type(n);
286 auto r0 = tensor_type(n);
287 auto r00 = tensor_complex_type(n);
288
289
290 auto one = value_type(1);
291 auto v = one;
292 for(auto& aa: a)
293 aa = v, v += one;
294
295 tensor_type b = (a+a) / value_type( 2 );
296 tensor_type r1 = ublas::real( (a+a) / value_type( 2 ) );
297 std::transform( b.begin(), b.end(), r0.begin(), [](auto const& l){ return std::real( l ); } );
298 BOOST_CHECK( r0 == r1 );
299
300 tensor_type r2 = ublas::imag( (a+a) / value_type( 2 ) );
301 std::transform( b.begin(), b.end(), r0.begin(), [](auto const& l){ return std::imag( l ); } );
302 BOOST_CHECK( r0 == r2 );
303
304 tensor_complex_type r3 = ublas::conj( (a+a) / value_type( 2 ) );
305 std::transform( b.begin(), b.end(), r00.begin(), [](auto const& l){ return std::conj( l ); } );
306 BOOST_CHECK( r00 == r3 );
307
308 }
309
310 for(auto const& n : extents) {
311
312
313
314
315 auto a = tensor_complex_type(n);
316
317 auto r00 = tensor_complex_type(n);
318 auto r0 = tensor_type(n);
319
320
321 auto one = complex_type(1,1);
322 auto v = one;
323 for(auto& aa: a)
324 aa = v, v = v + one;
325
326 tensor_complex_type b = (a+a) / complex_type( 2,2 );
327
328
329 tensor_type r1 = ublas::real( (a+a) / complex_type( 2,2 ) );
330 std::transform( b.begin(), b.end(), r0.begin(), [](auto const& l){ return std::real( l ); } );
331 BOOST_CHECK( r0 == r1 );
332
333 tensor_type r2 = ublas::imag( (a+a) / complex_type( 2,2 ) );
334 std::transform( b.begin(), b.end(), r0.begin(), [](auto const& l){ return std::imag( l ); } );
335 BOOST_CHECK( r0 == r2 );
336
337 tensor_complex_type r3 = ublas::conj( (a+a) / complex_type( 2,2 ) );
338 std::transform( b.begin(), b.end(), r00.begin(), [](auto const& l){ return std::conj( l ); } );
339 BOOST_CHECK( r00 == r3 );
340
341
342
343 }
344
345
346
347 }
348
349
350
351
352
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_outer_prod,value,test_types,fixture)353 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_outer_prod, value, test_types, fixture )
354 {
355 using namespace boost::numeric;
356 using value_type = typename value::first_type;
357 using layout_type = typename value::second_type;
358 using tensor_type = ublas::tensor<value_type,layout_type>;
359
360 for(auto const& n1 : extents) {
361 auto a = tensor_type(n1, value_type(2));
362 for(auto const& n2 : extents) {
363
364 auto b = tensor_type(n2, value_type(1));
365 auto c = ublas::outer_prod(a, b);
366
367 for(auto const& cc : c)
368 BOOST_CHECK_EQUAL( cc , a[0]*b[0] );
369 }
370 }
371 }
372
373
374
375 template<class V>
init(std::vector<V> & a)376 void init(std::vector<V>& a)
377 {
378 auto v = V(1);
379 for(auto i = 0u; i < a.size(); ++i, ++v){
380 a[i] = v;
381 }
382 }
383
384 template<class V>
init(std::vector<std::complex<V>> & a)385 void init(std::vector<std::complex<V>>& a)
386 {
387 auto v = std::complex<V>(1,1);
388 for(auto i = 0u; i < a.size(); ++i){
389 a[i] = v;
390 v.real(v.real()+1);
391 v.imag(v.imag()+1);
392 }
393 }
394
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_trans,value,test_types,fixture)395 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_trans, value, test_types, fixture )
396 {
397 using namespace boost::numeric;
398 using value_type = typename value::first_type;
399 using layout_type = typename value::second_type;
400 using tensor_type = ublas::tensor<value_type,layout_type>;
401
402 auto fak = [](auto const& p){
403 auto f = 1ul;
404 for(auto i = 1u; i <= p; ++i)
405 f *= i;
406 return f;
407 };
408
409 auto inverse = [](auto const& pi){
410 auto pi_inv = pi;
411 for(auto j = 0u; j < pi.size(); ++j)
412 pi_inv[pi[j]-1] = j+1;
413 return pi_inv;
414 };
415
416 for(auto const& n : extents)
417 {
418 auto const p = n.size();
419 auto const s = n.product();
420 auto aref = tensor_type(n);
421 auto v = value_type{};
422 for(auto i = 0u; i < s; ++i, v+=1)
423 aref[i] = v;
424 auto a = aref;
425
426
427 auto pi = std::vector<std::size_t>(p);
428 std::iota(pi.begin(), pi.end(), 1);
429 a = ublas::trans( a, pi );
430 BOOST_CHECK( a == aref );
431
432
433 auto const pfak = fak(p);
434 auto i = 0u;
435 for(; i < pfak-1; ++i) {
436 std::next_permutation(pi.begin(), pi.end());
437 a = ublas::trans( a, pi );
438 }
439 std::next_permutation(pi.begin(), pi.end());
440 for(; i > 0; --i) {
441 std::prev_permutation(pi.begin(), pi.end());
442 auto pi_inv = inverse(pi);
443 a = ublas::trans( a, pi_inv );
444 }
445
446 BOOST_CHECK( a == aref );
447
448 }
449 }
450
451
452 BOOST_AUTO_TEST_SUITE_END()
453
454