123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453 |
- // Copyright (c) 2018-2019 Cem Bassoy
- //
- // Distributed under the Boost Software License, Version 1.0. (See
- // accompanying file LICENSE_1_0.txt or copy at
- // http://www.boost.org/LICENSE_1_0.txt)
- //
- // The authors gratefully acknowledge the support of
- // Fraunhofer and Google in producing this work
- // which started as a Google Summer of Code project.
- //
- // And we acknowledge the support from all contributors.
- #include <iostream>
- #include <algorithm>
- #include <boost/numeric/ublas/tensor.hpp>
- #include <boost/numeric/ublas/matrix.hpp>
- #include <boost/numeric/ublas/vector.hpp>
- #include <boost/test/unit_test.hpp>
- #include "utility.hpp"
- BOOST_AUTO_TEST_SUITE ( test_tensor_functions, * boost::unit_test::depends_on("test_tensor_contraction") )
- using test_types = zip<int,long,float,double,std::complex<float>>::with_t<boost::numeric::ublas::first_order, boost::numeric::ublas::last_order>;
- //using test_types = zip<int>::with_t<boost::numeric::ublas::first_order>;
- struct fixture
- {
- using extents_type = boost::numeric::ublas::shape;
- fixture()
- : extents {
- extents_type{1,1}, // 1
- extents_type{1,2}, // 2
- extents_type{2,1}, // 3
- extents_type{2,3}, // 4
- extents_type{2,3,1}, // 5
- extents_type{4,1,3}, // 6
- extents_type{1,2,3}, // 7
- extents_type{4,2,3}, // 8
- extents_type{4,2,3,5}} // 9
- {
- }
- std::vector<extents_type> extents;
- };
- BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_vector, value, test_types, fixture )
- {
- using namespace boost::numeric;
- using value_type = typename value::first_type;
- using layout_type = typename value::second_type;
- using tensor_type = ublas::tensor<value_type,layout_type>;
- using vector_type = typename tensor_type::vector_type;
- for(auto const& n : extents){
- auto a = tensor_type(n, value_type{2});
- for(auto m = 0u; m < n.size(); ++m){
- auto b = vector_type (n[m], value_type{1} );
- auto c = ublas::prod(a, b, m+1);
- for(auto i = 0u; i < c.size(); ++i)
- BOOST_CHECK_EQUAL( c[i] , value_type(n[m]) * a[i] );
- }
- }
- }
- BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_matrix, value, test_types, fixture )
- {
- using namespace boost::numeric;
- using value_type = typename value::first_type;
- using layout_type = typename value::second_type;
- using tensor_type = ublas::tensor<value_type,layout_type>;
- using matrix_type = typename tensor_type::matrix_type;
- for(auto const& n : extents) {
- auto a = tensor_type(n, value_type{2});
- for(auto m = 0u; m < n.size(); ++m){
- auto b = matrix_type ( n[m], n[m], value_type{1} );
- auto c = ublas::prod(a, b, m+1);
- for(auto i = 0u; i < c.size(); ++i)
- BOOST_CHECK_EQUAL( c[i] , value_type(n[m]) * a[i] );
- }
- }
- }
- BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_tensor_1, value, test_types, fixture )
- {
- using namespace boost::numeric;
- using value_type = typename value::first_type;
- using layout_type = typename value::second_type;
- using tensor_type = ublas::tensor<value_type,layout_type>;
- // left-hand and right-hand side have the
- // the same number of elements
- for(auto const& na : extents) {
- auto a = tensor_type( na, value_type{2} );
- auto b = tensor_type( na, value_type{3} );
- auto const pa = a.rank();
- // the number of contractions is changed.
- for( auto q = 0ul; q <= pa; ++q) { // pa
- auto phi = std::vector<std::size_t> ( q );
- std::iota(phi.begin(), phi.end(), 1ul);
- auto c = ublas::prod(a, b, phi);
- auto acc = value_type(1);
- for(auto i = 0ul; i < q; ++i)
- acc *= a.extents().at(phi.at(i)-1);
- for(auto i = 0ul; i < c.size(); ++i)
- BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
- }
- }
- }
- BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_prod_tensor_2, value, test_types, fixture )
- {
- using namespace boost::numeric;
- using value_type = typename value::first_type;
- using layout_type = typename value::second_type;
- using tensor_type = ublas::tensor<value_type,layout_type>;
- auto compute_factorial = [](auto const& p){
- auto f = 1ul;
- for(auto i = 1u; i <= p; ++i)
- f *= i;
- return f;
- };
- auto permute_extents = [](auto const& pi, auto const& na){
- auto nb = na;
- assert(pi.size() == na.size());
- for(auto j = 0u; j < pi.size(); ++j)
- nb[pi[j]-1] = na[j];
- return nb;
- };
- // left-hand and right-hand side have the
- // the same number of elements
- for(auto const& na : extents) {
- auto a = tensor_type( na, value_type{2} );
- auto const pa = a.rank();
- auto pi = std::vector<std::size_t>(pa);
- auto fac = compute_factorial(pa);
- std::iota( pi.begin(), pi.end(), 1 );
- for(auto f = 0ul; f < fac; ++f)
- {
- auto nb = permute_extents( pi, na );
- auto b = tensor_type( nb, value_type{3} );
- // the number of contractions is changed.
- for( auto q = 0ul; q <= pa; ++q) { // pa
- auto phia = std::vector<std::size_t> ( q ); // concatenation for a
- auto phib = std::vector<std::size_t> ( q ); // concatenation for b
- std::iota(phia.begin(), phia.end(), 1ul);
- std::transform( phia.begin(), phia.end(), phib.begin(),
- [&pi] ( std::size_t i ) { return pi.at(i-1); } );
- auto c = ublas::prod(a, b, phia, phib);
- auto acc = value_type(1);
- for(auto i = 0ul; i < q; ++i)
- acc *= a.extents().at(phia.at(i)-1);
- for(auto i = 0ul; i < c.size(); ++i)
- BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
- }
- std::next_permutation(pi.begin(), pi.end());
- }
- }
- }
- BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_inner_prod, value, test_types, fixture )
- {
- using namespace boost::numeric;
- using value_type = typename value::first_type;
- using layout_type = typename value::second_type;
- using tensor_type = ublas::tensor<value_type,layout_type>;
- for(auto const& n : extents) {
- auto a = tensor_type(n, value_type(2));
- auto b = tensor_type(n, value_type(1));
- auto c = ublas::inner_prod(a, b);
- auto r = std::inner_product(a.begin(),a.end(), b.begin(),value_type(0));
- BOOST_CHECK_EQUAL( c , r );
- }
- }
- BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_norm, value, test_types, fixture )
- {
- using namespace boost::numeric;
- using value_type = typename value::first_type;
- using layout_type = typename value::second_type;
- using tensor_type = ublas::tensor<value_type,layout_type>;
- for(auto const& n : extents) {
- auto a = tensor_type(n);
- auto one = value_type(1);
- auto v = one;
- for(auto& aa: a)
- aa = v, v += one;
- auto c = ublas::inner_prod(a, a);
- auto r = std::inner_product(a.begin(),a.end(), a.begin(),value_type(0));
- auto r2 = ublas::norm( (a+a) / 2 );
- BOOST_CHECK_EQUAL( c , r );
- BOOST_CHECK_EQUAL( std::sqrt( c ) , r2 );
- }
- }
- BOOST_FIXTURE_TEST_CASE( test_tensor_real_imag_conj, fixture )
- {
- using namespace boost::numeric;
- using value_type = float;
- using complex_type = std::complex<value_type>;
- using layout_type = ublas::first_order;
- using tensor_complex_type = ublas::tensor<complex_type,layout_type>;
- using tensor_type = ublas::tensor<value_type,layout_type>;
- for(auto const& n : extents) {
- auto a = tensor_type(n);
- auto r0 = tensor_type(n);
- auto r00 = tensor_complex_type(n);
- auto one = value_type(1);
- auto v = one;
- for(auto& aa: a)
- aa = v, v += one;
- tensor_type b = (a+a) / value_type( 2 );
- tensor_type r1 = ublas::real( (a+a) / value_type( 2 ) );
- std::transform( b.begin(), b.end(), r0.begin(), [](auto const& l){ return std::real( l ); } );
- BOOST_CHECK( r0 == r1 );
- tensor_type r2 = ublas::imag( (a+a) / value_type( 2 ) );
- std::transform( b.begin(), b.end(), r0.begin(), [](auto const& l){ return std::imag( l ); } );
- BOOST_CHECK( r0 == r2 );
- tensor_complex_type r3 = ublas::conj( (a+a) / value_type( 2 ) );
- std::transform( b.begin(), b.end(), r00.begin(), [](auto const& l){ return std::conj( l ); } );
- BOOST_CHECK( r00 == r3 );
- }
- for(auto const& n : extents) {
- auto a = tensor_complex_type(n);
- auto r00 = tensor_complex_type(n);
- auto r0 = tensor_type(n);
- auto one = complex_type(1,1);
- auto v = one;
- for(auto& aa: a)
- aa = v, v = v + one;
- tensor_complex_type b = (a+a) / complex_type( 2,2 );
- tensor_type r1 = ublas::real( (a+a) / complex_type( 2,2 ) );
- std::transform( b.begin(), b.end(), r0.begin(), [](auto const& l){ return std::real( l ); } );
- BOOST_CHECK( r0 == r1 );
- tensor_type r2 = ublas::imag( (a+a) / complex_type( 2,2 ) );
- std::transform( b.begin(), b.end(), r0.begin(), [](auto const& l){ return std::imag( l ); } );
- BOOST_CHECK( r0 == r2 );
- tensor_complex_type r3 = ublas::conj( (a+a) / complex_type( 2,2 ) );
- std::transform( b.begin(), b.end(), r00.begin(), [](auto const& l){ return std::conj( l ); } );
- BOOST_CHECK( r00 == r3 );
- }
- }
- BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_outer_prod, value, test_types, fixture )
- {
- using namespace boost::numeric;
- using value_type = typename value::first_type;
- using layout_type = typename value::second_type;
- using tensor_type = ublas::tensor<value_type,layout_type>;
- for(auto const& n1 : extents) {
- auto a = tensor_type(n1, value_type(2));
- for(auto const& n2 : extents) {
- auto b = tensor_type(n2, value_type(1));
- auto c = ublas::outer_prod(a, b);
- for(auto const& cc : c)
- BOOST_CHECK_EQUAL( cc , a[0]*b[0] );
- }
- }
- }
- template<class V>
- void init(std::vector<V>& a)
- {
- auto v = V(1);
- for(auto i = 0u; i < a.size(); ++i, ++v){
- a[i] = v;
- }
- }
- template<class V>
- void init(std::vector<std::complex<V>>& a)
- {
- auto v = std::complex<V>(1,1);
- for(auto i = 0u; i < a.size(); ++i){
- a[i] = v;
- v.real(v.real()+1);
- v.imag(v.imag()+1);
- }
- }
- BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_trans, value, test_types, fixture )
- {
- using namespace boost::numeric;
- using value_type = typename value::first_type;
- using layout_type = typename value::second_type;
- using tensor_type = ublas::tensor<value_type,layout_type>;
- auto fak = [](auto const& p){
- auto f = 1ul;
- for(auto i = 1u; i <= p; ++i)
- f *= i;
- return f;
- };
- auto inverse = [](auto const& pi){
- auto pi_inv = pi;
- for(auto j = 0u; j < pi.size(); ++j)
- pi_inv[pi[j]-1] = j+1;
- return pi_inv;
- };
- for(auto const& n : extents)
- {
- auto const p = n.size();
- auto const s = n.product();
- auto aref = tensor_type(n);
- auto v = value_type{};
- for(auto i = 0u; i < s; ++i, v+=1)
- aref[i] = v;
- auto a = aref;
- auto pi = std::vector<std::size_t>(p);
- std::iota(pi.begin(), pi.end(), 1);
- a = ublas::trans( a, pi );
- BOOST_CHECK( a == aref );
- auto const pfak = fak(p);
- auto i = 0u;
- for(; i < pfak-1; ++i) {
- std::next_permutation(pi.begin(), pi.end());
- a = ublas::trans( a, pi );
- }
- std::next_permutation(pi.begin(), pi.end());
- for(; i > 0; --i) {
- std::prev_permutation(pi.begin(), pi.end());
- auto pi_inv = inverse(pi);
- a = ublas::trans( a, pi_inv );
- }
- BOOST_CHECK( a == aref );
- }
- }
- BOOST_AUTO_TEST_SUITE_END()
|