diff options
Diffstat (limited to 'src/boost/libs/numeric/ublas/test/tensor/test_functions.cpp')
| -rw-r--r-- | src/boost/libs/numeric/ublas/test/tensor/test_functions.cpp | 453 |
1 files changed, 453 insertions, 0 deletions
diff --git a/src/boost/libs/numeric/ublas/test/tensor/test_functions.cpp b/src/boost/libs/numeric/ublas/test/tensor/test_functions.cpp new file mode 100644 index 000000000..64ecda396 --- /dev/null +++ b/src/boost/libs/numeric/ublas/test/tensor/test_functions.cpp @@ -0,0 +1,453 @@ +// 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() + |