ceph/src/boost/libs/numeric/ublas/test/tensor/test_functions.cpp

   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;
  35         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
  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
  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
 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
 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
 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
 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
 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
 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>
 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>
 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
 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