ublas/test/tensor/test_multiplication.cpp
2019-02-25 17:05:02 +01:00

492 lines
13 KiB
C++

// 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.
//
#include <iostream>
#include <algorithm>
#include <vector>
#include <boost/numeric/ublas/tensor/multiplication.hpp>
#include <boost/numeric/ublas/tensor/extents.hpp>
#include <boost/numeric/ublas/tensor/strides.hpp>
#include "utility.hpp"
#include <boost/test/unit_test.hpp>
BOOST_AUTO_TEST_SUITE (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{5,4}, // 5
extents_type{2,3,1}, // 6
extents_type{4,1,3}, // 7
extents_type{1,2,3}, // 8
extents_type{4,2,3}, // 9
extents_type{4,2,3,5}} // 10
{
}
std::vector<extents_type> extents;
};
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_mtv, value, test_types, fixture )
{
using namespace boost::numeric;
using value_type = typename value::first_type;
using layout_type = typename value::second_type;
using strides_type = ublas::strides<layout_type>;
using vector_type = std::vector<value_type>;
using extents_type = ublas::shape;
using extents_type_base = typename extents_type::base_type;
using size_type = typename extents_type_base::value_type;
for(auto const& na : extents) {
if(na.size() > 2)
continue;
auto a = vector_type(na.product(), value_type{2});
auto wa = strides_type(na);
for(auto m = 0u; m < na.size(); ++m){
auto nb = extents_type {na[m],1};
auto wb = strides_type (nb);
auto b = vector_type (nb.product(), value_type{1} );
auto nc_base = extents_type_base(std::max(na.size()-1, size_type{2}), 1);
for(auto i = 0u, j = 0u; i < na.size(); ++i)
if(i != m)
nc_base[j++] = na[i];
auto nc = extents_type (nc_base);
auto wc = strides_type (nc);
auto c = vector_type (nc.product(), value_type{0});
ublas::detail::recursive::mtv(
size_type(m),
c.data(), nc.data(), wc.data(),
a.data(), na.data(), wa.data(),
b.data());
for(auto i = 0u; i < c.size(); ++i)
BOOST_CHECK_EQUAL( c[i] , value_type(na[m]) * a[i] );
}
}
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_mtm, value, test_types, fixture )
{
using namespace boost::numeric;
using value_type = typename value::first_type;
using layout_type = typename value::second_type;
using strides_type = ublas::strides<layout_type>;
using vector_type = std::vector<value_type>;
using extents_type = ublas::shape;
// using extents_type_base = typename extents_type::base_type;
for(auto const& na : extents) {
if(na.size() != 2)
continue;
auto a = vector_type (na.product(), value_type{2});
auto wa = strides_type (na);
auto nb = extents_type {na[1],na[0]};
auto wb = strides_type (nb);
auto b = vector_type (nb.product(), value_type{1} );
auto nc = extents_type {na[0],nb[1]};
auto wc = strides_type (nc);
auto c = vector_type (nc.product());
ublas::detail::recursive::mtm(
c.data(), nc.data(), wc.data(),
a.data(), na.data(), wa.data(),
b.data(), nb.data(), wb.data());
for(auto i = 0u; i < c.size(); ++i)
BOOST_CHECK_EQUAL( c[i] , value_type(na[1]) * a[0] );
}
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttv, value, test_types, fixture )
{
using namespace boost::numeric;
using value_type = typename value::first_type;
using layout_type = typename value::second_type;
using strides_type = ublas::strides<layout_type>;
using vector_type = std::vector<value_type>;
using extents_type = ublas::shape;
using extents_type_base = typename extents_type::base_type;
using size_type = typename extents_type_base::value_type;
for(auto const& na : extents) {
auto a = vector_type(na.product(), value_type{2});
auto wa = strides_type(na);
for(auto m = 0u; m < na.size(); ++m){
auto b = vector_type (na[m], value_type{1} );
auto nb = extents_type {na[m],1};
auto wb = strides_type (nb);
auto nc_base = extents_type_base(std::max(na.size()-1, size_type(2)),1);
for(auto i = 0ul, j = 0ul; i < na.size(); ++i)
if(i != m)
nc_base[j++] = na[i];
auto nc = extents_type (nc_base);
auto wc = strides_type (nc);
auto c = vector_type (nc.product(), value_type{0});
ublas::ttv(size_type(m+1), na.size(),
c.data(), nc.data(), wc.data(),
a.data(), na.data(), wa.data(),
b.data(), nb.data(), wb.data());
for(auto i = 0u; i < c.size(); ++i)
BOOST_CHECK_EQUAL( c[i] , value_type(na[m]) * a[i] );
}
}
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttm, value, test_types, fixture )
{
using namespace boost::numeric;
using value_type = typename value::first_type;
using layout_type = typename value::second_type;
using strides_type = ublas::strides<layout_type>;
using vector_type = std::vector<value_type>;
using extents_type = ublas::shape;
using size_type = typename extents_type::value_type;
for(auto const& na : extents) {
auto a = vector_type(na.product(), value_type{2});
auto wa = strides_type(na);
for(auto m = 0u; m < na.size(); ++m){
auto nb = extents_type {na[m], na[m] };
auto b = vector_type (nb.product(), value_type{1} );
auto wb = strides_type (nb);
auto nc = na;
auto wc = strides_type (nc);
auto c = vector_type (nc.product(), value_type{0});
ublas::ttm(size_type(m+1), na.size(),
c.data(), nc.data(), wc.data(),
a.data(), na.data(), wa.data(),
b.data(), nb.data(), wb.data());
for(auto i = 0u; i < c.size(); ++i)
BOOST_CHECK_EQUAL( c[i] , value_type(na[m]) * a[i] );
}
}
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttt_permutation, value, test_types, fixture )
{
using namespace boost::numeric;
using value_type = typename value::first_type;
using layout_type = typename value::second_type;
using strides_type = ublas::strides<layout_type>;
using vector_type = std::vector<value_type>;
using extents_type = ublas::shape;
using size_type = typename strides_type::value_type;
auto compute_factorial = [](auto const& p){
auto f = 1ul;
for(auto i = 1u; i <= p; ++i)
f *= i;
return f;
};
auto compute_inverse_permutation = [](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;
};
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[j] = na[pi[j]-1];
return nb;
};
// left-hand and right-hand side have the
// the same number of elements
// computing the inner product with
// different permutation tuples for
// right-hand side
for(auto const& na : extents) {
auto wa = strides_type(na);
auto a = vector_type(na.product(), value_type{2});
auto pa = na.size();
auto pia = std::vector<size_type>(pa);
std::iota( pia.begin(), pia.end(), 1 );
auto pib = pia;
auto pib_inv = compute_inverse_permutation(pib);
auto f = compute_factorial(pa);
// for the number of possible permutations
// only permutation tuple pib is changed.
for(auto i = 0u; i < f; ++i) {
auto nb = permute_extents( pib, na );
auto wb = strides_type(nb);
auto b = vector_type(nb.product(), value_type{3});
auto pb = nb.size();
// the number of contractions is changed.
for( auto q = size_type(0); q <= pa; ++q) {
auto r = pa - q;
auto s = pb - q;
auto pc = r+s > 0 ? std::max(r+s,size_type(2)) : size_type(2);
auto nc_base = std::vector<size_type>( pc , 1 );
for(auto i = 0u; i < r; ++i)
nc_base[ i ] = na[ pia[i]-1 ];
for(auto i = 0u; i < s; ++i)
nc_base[ r + i ] = nb[ pib_inv[i]-1 ];
auto nc = extents_type ( nc_base );
auto wc = strides_type ( nc );
auto c = vector_type ( nc.product(), value_type(0) );
ublas::ttt(pa,pb,q,
pia.data(), pib_inv.data(),
c.data(), nc.data(), wc.data(),
a.data(), na.data(), wa.data(),
b.data(), nb.data(), wb.data());
auto acc = value_type(1);
for(auto i = r; i < pa; ++i)
acc *= value_type(na[pia[i]-1]);
for(auto i = 0ul; i < c.size(); ++i)
BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
}
std::next_permutation(pib.begin(), pib.end());
pib_inv = compute_inverse_permutation(pib);
}
}
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttt, value, test_types, fixture )
{
using namespace boost::numeric;
using value_type = typename value::first_type;
using layout_type = typename value::second_type;
using strides_type = ublas::strides<layout_type>;
using vector_type = std::vector<value_type>;
using extents_type = ublas::shape;
using size_type = typename strides_type::value_type;
// left-hand and right-hand side have the
// the same number of elements
// computing the inner product with
// different permutation tuples for
// right-hand side
for(auto const& na : extents) {
auto wa = strides_type(na);
auto a = vector_type(na.product(), value_type{2});
auto pa = na.size();
auto nb = na;
auto wb = strides_type(nb);
auto b = vector_type(nb.product(), value_type{3});
auto pb = nb.size();
// std::cout << "na = ";
// std::copy(na.begin(), na.end(), std::ostream_iterator<size_type>(std::cout, " "));
// std::cout << std::endl;
// std::cout << "nb = ";
// std::copy(nb.begin(), nb.end(), std::ostream_iterator<size_type>(std::cout, " "));
// std::cout << std::endl;
// the number of contractions is changed.
for( auto q = size_type(0); q <= pa; ++q) { // pa
auto r = pa - q;
auto s = pb - q;
auto pc = r+s > 0 ? std::max(r+s, size_type(2)) : size_type(2);
auto nc_base = std::vector<size_type>( pc , 1 );
for(auto i = 0u; i < r; ++i)
nc_base[ i ] = na[ i ];
for(auto i = 0u; i < s; ++i)
nc_base[ r + i ] = nb[ i ];
auto nc = extents_type ( nc_base );
auto wc = strides_type ( nc );
auto c = vector_type ( nc.product(), value_type{0} );
// std::cout << "nc = ";
// std::copy(nc.begin(), nc.end(), std::ostream_iterator<size_type>(std::cout, " "));
// std::cout << std::endl;
ublas::ttt(pa,pb,q,
c.data(), nc.data(), wc.data(),
a.data(), na.data(), wa.data(),
b.data(), nb.data(), wb.data());
auto acc = value_type(1);
for(auto i = r; i < pa; ++i)
acc *= value_type(na[i]);
for(auto i = 0u; i < c.size(); ++i)
BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
}
}
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_inner, value, test_types, fixture )
{
using namespace boost::numeric;
using value_type = typename value::first_type;
using layout_type = typename value::second_type;
using strides_type = ublas::strides<layout_type>;
using vector_type = std::vector<value_type>;
for(auto const& n : extents) {
auto a = vector_type(n.product(), value_type{2});
auto b = vector_type(n.product(), value_type{3});
auto w = strides_type(n);
auto c = ublas::inner(n.size(), n.data(), a.data(), w.data(), b.data(), w.data(), value_type(0));
auto cref = std::inner_product(a.begin(), a.end(), b.begin(), value_type(0));
BOOST_CHECK_EQUAL( c , cref );
}
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_outer, value, test_types, fixture )
{
using namespace boost::numeric;
using value_type = typename value::first_type;
using layout_type = typename value::second_type;
using extents_type = ublas::shape;
using strides_type = ublas::strides<layout_type>;
using vector_type = std::vector<value_type>;
for(auto const& na : extents) {
auto a = vector_type(na.product(), value_type{2});
auto wa = strides_type(na);
for(auto const& nb : extents) {
auto b = vector_type(nb.product(), value_type{3});
auto wb = strides_type(nb);
auto c = vector_type(nb.product()*na.product());
auto nc = typename extents_type::base_type(na.size()+nb.size());
for(auto i = 0u; i < na.size(); ++i)
nc[i] = na[i];
for(auto i = 0u; i < nb.size(); ++i)
nc[i+na.size()] = nb[i];
auto wc = strides_type(extents_type(nc));
ublas::outer(c.data(), nc.size(), nc.data(), wc.data(),
a.data(), na.size(), na.data(), wa.data(),
b.data(), nb.size(), nb.data(), wb.data());
for(auto const& cc : c)
BOOST_CHECK_EQUAL( cc , a[0]*b[0] );
}
}
}
BOOST_AUTO_TEST_SUITE_END()