ublas/tensor/functions.hpp

0001 //
0002 //  Copyright (c) 2018-2019, Cem Bassoy, cem.bassoy@gmail.com
0003 //
0004 //  Distributed under the Boost Software License, Version 1.0. (See
0005 //  accompanying file LICENSE_1_0.txt or copy at
0006 //  http://www.boost.org/LICENSE_1_0.txt)
0007 //
0008 //  The authors gratefully acknowledge the support of
0009 //  Fraunhofer IOSB, Ettlingen, Germany
0010 //
0011
0012
0013 #ifndef BOOST_UBLAS_TENSOR_FUNCTIONS_HPP
0014 #define BOOST_UBLAS_TENSOR_FUNCTIONS_HPP
0015
0016
0017 #include <stdexcept>
0018 #include <vector>
0019 #include <algorithm>
0020 #include <numeric>
0021
0022
0023 #include "multiplication.hpp"
0024 #include "algorithms.hpp"
0025 #include "expression.hpp"
0026 #include "expression_evaluation.hpp"
0027 #include "storage_traits.hpp"
0028
0029 namespace boost {
0030 namespace numeric {
0031 namespace ublas {
0032
0033 template<class Value, class Format, class Allocator>
0034 class tensor;
0035
0036 template<class Value, class Format, class Allocator>
0037 class matrix;
0038
0039 template<class Value, class Allocator>
0040 class vector;
0041
0042
0043
0044
0045 /** @brief Computes the m-mode tensor-times-vector product
0046  *
0047  * Implements C[i1,...,im-1,im+1,...,ip] = A[i1,i2,...,ip] * b[im]
0048  *
0049  * @note calls ublas::ttv
0050  *
0051  * @param[in] m contraction dimension with 1 <= m <= p
0052  * @param[in] a tensor object A with order p
0053  * @param[in] b vector object B
0054  *
0055  * @returns tensor object C with order p-1, the same storage format and allocator type as A
0056 */
0057 template<class V, class F, class A1, class A2>
0058 auto prod(tensor<V,F,A1> const& a, vector<V,A2> const& b, const std::size_t m)
0059 {
0060
0061     using tensor_type  = tensor<V,F,A1>;
0062     using extents_type = typename tensor_type::extents_type;
0063     using ebase_type   = typename extents_type::base_type;
0064     using value_type   = typename tensor_type::value_type;
0065     using size_type = typename extents_type::value_type;
0066
0067     auto const p = std::size_t(a.rank());
0068
0069     if( m == 0)
0070         throw std::length_error("error in boost::numeric::ublas::prod(ttv): contraction mode must be greater than zero.");
0071
0072     if( p < m )
0073         throw std::length_error("error in boost::numeric::ublas::prod(ttv): rank of tensor must be greater than or equal to the modus.");
0074
0075     if( p == 0)
0076         throw std::length_error("error in boost::numeric::ublas::prod(ttv): rank of tensor must be greater than zero.");
0077
0078     if( a.empty() )
0079         throw std::length_error("error in boost::numeric::ublas::prod(ttv): first argument tensor should not be empty.");
0080
0081     if( b.size() == 0)
0082         throw std::length_error("error in boost::numeric::ublas::prod(ttv): second argument vector should not be empty.");
0083
0084
0085     auto nc = ebase_type(std::max(p-1, size_type(2)) , size_type(1));
0086     auto nb = ebase_type{b.size(),1};
0087
0088
0089     for(auto i = 0u, j = 0u; i < p; ++i)
0090         if(i != m-1)
0091             nc[j++] = a.extents().at(i);
0092
0093     auto c = tensor_type(extents_type(nc),value_type{});
0094
0095     auto bb = &(b(0));
0096
0097     ttv(m, p,
0098         c.data(), c.extents().data(), c.strides().data(),
0099         a.data(), a.extents().data(), a.strides().data(),
0100         bb, nb.data(), nb.data());
0101
0102
0103     return c;
0104 }
0105
0106
0107
0108 /** @brief Computes the m-mode tensor-times-matrix product
0109  *
0110  * Implements C[i1,...,im-1,j,im+1,...,ip] = A[i1,i2,...,ip] * B[j,im]
0111  *
0112  * @note calls ublas::ttm
0113  *
0114  * @param[in] a tensor object A with order p
0115  * @param[in] b vector object B
0116  * @param[in] m contraction dimension with 1 <= m <= p
0117  *
0118  * @returns tensor object C with order p, the same storage format and allocator type as A
0119 */
0120 template<class V, class F, class A1, class A2>
0121 auto prod(tensor<V,F,A1> const& a, matrix<V,F,A2> const& b, const std::size_t m)
0122 {
0123
0124     using tensor_type  = tensor<V,F,A1>;
0125     using extents_type = typename tensor_type::extents_type;
0126     using strides_type = typename tensor_type::strides_type;
0127     using value_type   = typename tensor_type::value_type;
0128
0129
0130     auto const p = a.rank();
0131
0132     if( m == 0)
0133         throw std::length_error("error in boost::numeric::ublas::prod(ttm): contraction mode must be greater than zero.");
0134
0135     if( p < m || m > a.extents().size())
0136         throw std::length_error("error in boost::numeric::ublas::prod(ttm): rank of the tensor must be greater equal the modus.");
0137
0138     if( p == 0)
0139         throw std::length_error("error in boost::numeric::ublas::prod(ttm): rank of the tensor must be greater than zero.");
0140
0141     if( a.empty() )
0142         throw std::length_error("error in boost::numeric::ublas::prod(ttm): first argument tensor should not be empty.");
0143
0144     if( b.size1()*b.size2() == 0)
0145         throw std::length_error("error in boost::numeric::ublas::prod(ttm): second argument matrix should not be empty.");
0146
0147
0148     auto nc = a.extents().base();
0149     auto nb = extents_type {b.size1(),b.size2()};
0150     auto wb = strides_type (nb);
0151
0152     nc[m-1] = nb[0];
0153
0154     auto c = tensor_type(extents_type(nc),value_type{});
0155
0156     auto bb = &(b(0,0));
0157
0158     ttm(m, p,
0159         c.data(), c.extents().data(), c.strides().data(),
0160         a.data(), a.extents().data(), a.strides().data(),
0161         bb, nb.data(), wb.data());
0162
0163
0164     return c;
0165 }
0166
0167
0168
0169
0170 /** @brief Computes the q-mode tensor-times-tensor product
0171  *
0172  * Implements C[i1,...,ir,j1,...,js] = sum( A[i1,...,ir+q] * B[j1,...,js+q]  )
0173  *
0174  * @note calls ublas::ttt
0175  *
0176  * na[phia[x]] = nb[phib[x]] for 1 <= x <= q
0177  *
0178  * @param[in]    phia one-based permutation tuple of length q for the first input tensor a
0179  * @param[in]    phib one-based permutation tuple of length q for the second input tensor b
0180  * @param[in]  a  left-hand side tensor with order r+q
0181  * @param[in]  b  right-hand side tensor with order s+q
0182  * @result     tensor with order r+s
0183 */
0184 template<class V, class F, class A1, class A2>
0185 auto prod(tensor<V,F,A1> const& a, tensor<V,F,A2> const& b,
0186           std::vector<std::size_t> const& phia, std::vector<std::size_t> const& phib)
0187 {
0188
0189     using tensor_type  = tensor<V,F,A1>;
0190     using extents_type = typename tensor_type::extents_type;
0191     using value_type   = typename tensor_type::value_type;
0192     using size_type = typename extents_type::value_type;
0193
0194     auto const pa = a.rank();
0195     auto const pb = b.rank();
0196
0197     auto const q  = size_type(phia.size());
0198
0199     if(pa == 0ul)
0200         throw std::runtime_error("error in ublas::prod: order of left-hand side tensor must be greater than 0.");
0201     if(pb == 0ul)
0202         throw std::runtime_error("error in ublas::prod: order of right-hand side tensor must be greater than 0.");
0203     if(pa < q)
0204         throw std::runtime_error("error in ublas::prod: number of contraction dimensions cannot be greater than the order of the left-hand side tensor.");
0205     if(pb < q)
0206         throw std::runtime_error("error in ublas::prod: number of contraction dimensions cannot be greater than the order of the right-hand side tensor.");
0207
0208     if(q != phib.size())
0209         throw std::runtime_error("error in ublas::prod: permutation tuples must have the same length.");
0210
0211     if(pa < phia.size())
0212         throw std::runtime_error("error in ublas::prod: permutation tuple for the left-hand side tensor cannot be greater than the corresponding order.");
0213     if(pb < phib.size())
0214         throw std::runtime_error("error in ublas::prod: permutation tuple for the right-hand side tensor cannot be greater than the corresponding order.");
0215
0216
0217     auto const& na = a.extents();
0218     auto const& nb = b.extents();
0219
0220     for(auto i = 0ul; i < q; ++i)
0221         if( na.at(phia.at(i)-1) != nb.at(phib.at(i)-1))
0222             throw std::runtime_error("error in ublas::prod: permutations of the extents are not correct.");
0223
0224     auto const r = pa - q;
0225     auto const s = pb - q;
0226
0227
0228     std::vector<std::size_t> phia1(pa), phib1(pb);
0229     std::iota(phia1.begin(), phia1.end(), 1ul);
0230     std::iota(phib1.begin(), phib1.end(), 1ul);
0231
0232     std::vector<std::size_t> nc( std::max ( r+s , size_type(2) ), size_type(1) );
0233
0234     for(auto i = 0ul; i < phia.size(); ++i)
0235         * std::remove(phia1.begin(), phia1.end(), phia.at(i)) = phia.at(i);
0236
0237     //phia1.erase( std::remove(phia1.begin(), phia1.end(), phia.at(i)),  phia1.end() )  ;
0238
0239     assert(phia1.size() == pa);
0240
0241     for(auto i = 0ul; i < r; ++i)
0242         nc[ i ] = na[ phia1[ i  ] - 1  ];
0243
0244     for(auto i = 0ul; i < phib.size(); ++i)
0245         * std::remove(phib1.begin(), phib1.end(), phib.at(i))  = phib.at(i) ;
0246     //phib1.erase( std::remove(phib1.begin(), phib1.end(), phia.at(i)), phib1.end() )  ;
0247
0248     assert(phib1.size() == pb);
0249
0250     for(auto i = 0ul; i < s; ++i)
0251         nc[ r + i ] = nb[ phib1[ i  ] - 1  ];
0252
0253     //  std::copy( phib.begin(), phib.end(), phib1.end()  );
0254
0255     assert(  phia1.size() == pa  );
0256     assert(  phib1.size() == pb  );
0257
0258     auto c = tensor_type(extents_type(nc), value_type{});
0259
0260     ttt(pa, pb, q,
0261         phia1.data(), phib1.data(),
0262         c.data(), c.extents().data(), c.strides().data(),
0263         a.data(), a.extents().data(), a.strides().data(),
0264         b.data(), b.extents().data(), b.strides().data());
0265
0266     return c;
0267 }
0268
0269 //template<class V, class F, class A1, class A2, std::size_t N, std::size_t M>
0270 //auto operator*( tensor_index<V,F,A1,N> const& lhs, tensor_index<V,F,A2,M> const& rhs)
0271
0272
0273
0274
0275 /** @brief Computes the q-mode tensor-times-tensor product
0276  *
0277  * Implements C[i1,...,ir,j1,...,js] = sum( A[i1,...,ir+q] * B[j1,...,js+q]  )
0278  *
0279  * @note calls ublas::ttt
0280  *
0281  * na[phi[x]] = nb[phi[x]] for 1 <= x <= q
0282  *
0283  * @param[in]    phi one-based permutation tuple of length q for bot input tensors
0284  * @param[in]  a  left-hand side tensor with order r+q
0285  * @param[in]  b  right-hand side tensor with order s+q
0286  * @result     tensor with order r+s
0287 */
0288 template<class V, class F, class A1, class A2>
0289 auto prod(tensor<V,F,A1> const& a, tensor<V,F,A2> const& b,
0290           std::vector<std::size_t> const& phi)
0291 {
0292     return prod(a, b, phi, phi);
0293 }
0294
0295
0296 /** @brief Computes the inner product of two tensors
0297  *
0298  * Implements c = sum(A[i1,i2,...,ip] * B[i1,i2,...,jp])
0299  *
0300  * @note calls inner function
0301  *
0302  * @param[in] a tensor object A
0303  * @param[in] b tensor object B
0304  *
0305  * @returns a value type.
0306 */
0307 template<class V, class F, class A1, class A2>
0308 auto inner_prod(tensor<V,F,A1> const& a, tensor<V,F,A2> const& b)
0309 {
0310     using value_type   = typename tensor<V,F,A1>::value_type;
0311
0312     if( a.rank() != b.rank() )
0313         throw std::length_error("error in boost::numeric::ublas::inner_prod: Rank of both tensors must be the same.");
0314
0315     if( a.empty() || b.empty())
0316         throw std::length_error("error in boost::numeric::ublas::inner_prod: Tensors should not be empty.");
0317
0318     if( a.extents() != b.extents())
0319         throw std::length_error("error in boost::numeric::ublas::inner_prod: Tensor extents should be the same.");
0320
0321     return inner(a.rank(), a.extents().data(),
0322                  a.data(), a.strides().data(),
0323                  b.data(), b.strides().data(), value_type{0});
0324 }
0325
0326 /** @brief Computes the outer product of two tensors
0327  *
0328  * Implements C[i1,...,ip,j1,...,jq] = A[i1,i2,...,ip] * B[j1,j2,...,jq]
0329  *
0330  * @note calls outer function
0331  *
0332  * @param[in] a tensor object A
0333  * @param[in] b tensor object B
0334  *
0335  * @returns tensor object C with the same storage format F and allocator type A1
0336 */
0337 template<class V, class F, class A1, class A2>
0338 auto outer_prod(tensor<V,F,A1> const& a, tensor<V,F,A2> const& b)
0339 {
0340     using tensor_type  = tensor<V,F,A1>;
0341     using extents_type = typename tensor_type::extents_type;
0342
0343     if( a.empty() || b.empty() )
0344         throw std::runtime_error("error in boost::numeric::ublas::outer_prod: tensors should not be empty.");
0345
0346     auto nc = typename extents_type::base_type(a.rank() + b.rank());
0347     for(auto i = 0u; i < a.rank(); ++i)
0348         nc.at(i) = a.extents().at(i);
0349
0350     for(auto i = 0u; i < b.rank(); ++i)
0351         nc.at(a.rank()+i) = b.extents().at(i);
0352
0353     auto c = tensor_type(extents_type(nc));
0354
0355     outer(c.data(), c.rank(), c.extents().data(), c.strides().data(),
0356           a.data(), a.rank(), a.extents().data(), a.strides().data(),
0357           b.data(), b.rank(), b.extents().data(), b.strides().data());
0358
0359     return c;
0360 }
0361
0362
0363
0364 /** @brief Transposes a tensor according to a permutation tuple
0365  *
0366  * Implements C[tau[i1],tau[i2]...,tau[ip]] = A[i1,i2,...,ip]
0367  *
0368  * @note calls trans function
0369  *
0370  * @param[in] a    tensor object of rank p
0371  * @param[in] tau  one-based permutation tuple of length p
0372  * @returns        a transposed tensor object with the same storage format F and allocator type A
0373 */
0374 template<class V, class F, class A>
0375 auto trans(tensor<V,F,A> const& a, std::vector<std::size_t> const& tau)
0376 {
0377     using tensor_type  = tensor<V,F,A>;
0378     using extents_type = typename tensor_type::extents_type;
0379     //  using strides_type = typename tensor_type::strides_type;
0380
0381     if( a.empty() )
0382         return tensor<V,F,A>{};
0383
0384     auto const   p = a.rank();
0385     auto const& na = a.extents();
0386
0387     auto nc = typename extents_type::base_type (p);
0388     for(auto i = 0u; i < p; ++i)
0389         nc.at(tau.at(i)-1) = na.at(i);
0390
0391     //  auto wc = strides_type(extents_type(nc));
0392
0393     auto c = tensor_type(extents_type(nc));
0394
0395
0396     trans( a.rank(), a.extents().data(), tau.data(),
0397            c.data(), c.strides().data(),
0398            a.data(), a.strides().data());
0399
0400     //  auto wc_pi = typename strides_type::base_type (p);
0401     //  for(auto i = 0u; i < p; ++i)
0402     //      wc_pi.at(tau.at(i)-1) = c.strides().at(i);
0403
0404
0405     //copy(a.rank(),
0406     //       a.extents().data(),
0407     //       c.data(), wc_pi.data(),
0408     //       a.data(), a.strides().data() );
0409
0410     return c;
0411 }
0412
0413 /** @brief Computes the frobenius norm of a tensor expression
0414  *
0415  * @note evaluates the tensor expression and calls the accumulate function
0416  *
0417  *
0418  * Implements the two-norm with
0419  * k = sqrt( sum_(i1,...,ip) A(i1,...,ip)^2 )
0420  *
0421  * @param[in] a    tensor object of rank p
0422  * @returns        the frobenius norm of the tensor
0423 */
0424 //template<class V, class F, class A>
0425 //auto norm(tensor<V,F,A> const& a)
0426 template<class T, class D>
0427 auto norm(detail::tensor_expression<T,D> const& expr)
0428 {
0429
0430     using tensor_type = typename detail::tensor_expression<T,D>::tensor_type;
0431     using value_type = typename tensor_type::value_type;
0432
0433     auto a = tensor_type( expr );
0434
0435     if( a.empty() )
0436         throw std::runtime_error("error in boost::numeric::ublas::norm: tensors should not be empty.");
0437
0438     return std::sqrt( accumulate( a.order(), a.extents().data(), a.data(), a.strides().data(), value_type{},
0439                                   [](auto const& l, auto const& r){ return l + r*r; }  ) ) ;
0440 }
0441
0442
0443
0444 /** @brief Extract the real component of tensor elements within a tensor expression
0445  *
0446  * @param[in] lhs tensor expression
0447  * @returns   unary tensor expression
0448 */
0449 template<class T, class D>
0450 auto real(detail::tensor_expression<T,D> const& expr) {
0451     return detail::make_unary_tensor_expression<T> (expr(), [] (auto const& l) { return std::real( l ); } );
0452 }
0453
0454 /** @brief Extract the real component of tensor elements within a tensor expression
0455  *
0456  * @param[in] lhs tensor expression
0457  * @returns   unary tensor expression
0458 */
0459 template<class V, class F, class A, class D>
0460 auto real(detail::tensor_expression<tensor<std::complex<V>,F,A>,D> const& expr)
0461 {
0462     using tensor_complex_type = tensor<std::complex<V>,F,A>;
0463     using tensor_type = tensor<V,F,typename storage_traits<A>::template rebind<V>>;
0464
0465     if( detail::retrieve_extents( expr  ).empty() )
0466         throw std::runtime_error("error in boost::numeric::ublas::real: tensors should not be empty.");
0467
0468     auto a = tensor_complex_type( expr );
0469     auto c = tensor_type( a.extents() );
0470
0471     std::transform( a.begin(), a.end(),  c.begin(), [](auto const& l){ return std::real(l) ; }  );
0472
0473     return c;
0474 }
0475
0476
0477 /** @brief Extract the imaginary component of tensor elements within a tensor expression
0478  *
0479  * @param[in] lhs tensor expression
0480  * @returns   unary tensor expression
0481 */
0482 template<class T, class D>
0483 auto imag(detail::tensor_expression<T,D> const& lhs) {
0484     return detail::make_unary_tensor_expression<T> (lhs(), [] (auto const& l) { return std::imag( l ); } );
0485 }
0486
0487
0488 /** @brief Extract the imag component of tensor elements within a tensor expression
0489  *
0490  * @param[in] lhs tensor expression
0491  * @returns   unary tensor expression
0492 */
0493 template<class V, class A, class F, class D>
0494 auto imag(detail::tensor_expression<tensor<std::complex<V>,F,A>,D> const& expr)
0495 {
0496     using tensor_complex_type = tensor<std::complex<V>,F,A>;
0497     using tensor_type = tensor<V,F,typename storage_traits<A>::template rebind<V>>;
0498
0499     if( detail::retrieve_extents( expr  ).empty() )
0500         throw std::runtime_error("error in boost::numeric::ublas::real: tensors should not be empty.");
0501
0502     auto a = tensor_complex_type( expr );
0503     auto c = tensor_type( a.extents() );
0504
0505     std::transform( a.begin(), a.end(),  c.begin(), [](auto const& l){ return std::imag(l) ; }  );
0506
0507     return c;
0508 }
0509
0510 /** @brief Computes the complex conjugate component of tensor elements within a tensor expression
0511  *
0512  * @param[in] expr tensor expression
0513  * @returns   complex tensor
0514 */
0515 template<class T, class D>
0516 auto conj(detail::tensor_expression<T,D> const& expr)
0517 {
0518     using tensor_type = T;
0519     using value_type = typename tensor_type::value_type;
0520     using layout_type = typename tensor_type::layout_type;
0521     using array_type = typename tensor_type::array_type;
0522
0523     using new_value_type = std::complex<value_type>;
0524     using new_array_type = typename storage_traits<array_type>::template rebind<new_value_type>;
0525
0526     using tensor_complex_type = tensor<new_value_type,layout_type, new_array_type>;
0527
0528     if( detail::retrieve_extents( expr  ).empty() )
0529         throw std::runtime_error("error in boost::numeric::ublas::conj: tensors should not be empty.");
0530
0531     auto a = tensor_type( expr );
0532     auto c = tensor_complex_type( a.extents() );
0533
0534     std::transform( a.begin(), a.end(),  c.begin(), [](auto const& l){ return std::conj(l) ; }  );
0535
0536     return c;
0537 }
0538
0539
0540 /** @brief Computes the complex conjugate component of tensor elements within a tensor expression
0541  *
0542  * @param[in] lhs tensor expression
0543  * @returns   unary tensor expression
0544 */
0545 template<class V, class A, class F, class D>
0546 auto conj(detail::tensor_expression<tensor<std::complex<V>,F,A>,D> const& expr)
0547 {
0548     return detail::make_unary_tensor_expression<tensor<std::complex<V>,F,A>> (expr(), [] (auto const& l) { return std::conj( l ); } );
0549 }
0550
0551
0552
0553 }
0554 }
0555 }
0556
0557
0558 #endif