tensor.h
説明を見る。
1 //
2 // Copyright (c) 2003-2011, MIST Project, Intelligent Media Integration COE, Nagoya University
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification,
6 // are permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the Nagoya University nor the names of its contributors
16 // may be used to endorse or promote products derived from this software
17 // without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28 
29 
30 
37 #ifndef __INCLUDE_MIST_TENSOR__
38 #define __INCLUDE_MIST_TENSOR__
39 
40 #ifndef __INCLUDE_MIST_H__
41 #include "mist.h"
42 #endif
43 
44 #ifndef __INCLUDE_MIST_MATRIX__
45 #include "matrix.h"
46 #endif
47 
48 #ifndef __INCLUDE_NUMERIC__
49 #include "numeric.h"
50 #endif
51 
52 #include <numeric>
53 #include <functional>
54 
55 
56 // beginning of mist namespace
57 _MIST_BEGIN
58 
59 
60 template< int M, typename V, typename A >
61 class tensor;
62 
63 namespace __tensor__
64 {
65  template< int M, typename V, typename A >
66  class _base_
67  {
68  protected:
69  typedef V value_type;
70  typedef A allocator;
71  typedef typename array< value_type, allocator >::size_type size_type;
72  typedef matrix< V, A > matrix_type;
73 
74  size_type _mode;
75  array< size_type > _ranks;
76  array< value_type, allocator > _data;
77 
78  template< typename MV, typename MA >
79  void _fold( const size_type &m, const matrix< MV, MA > &mat )
80  {
81  array< size_type > i( _mode ), d( _mode );
82  for( size_type ii = 0 ; ii < _mode ; d[ ii ] = _diff( ii ), ii ++ );
83  for( typename matrix_type::size_type c = 0 ; c < mat.cols( ) ; c ++ )
84  {
85  for( typename matrix_type::size_type r = i[ m - 1 ] = 0 ; r < mat.rows( ) ; r ++, i[ m - 1 ] ++ )
86  {
87  operator [ ]( ::std::inner_product( &i[ 0 ], &i[ 0 ] + i.size( ), &d[ 0 ], 0 ) ) = mat( r, c );
88  }
89  _inc( m % _mode, m - 1, i );
90  }
91  }
92 
93  private:
94  size_type _diff( const size_type &i ) const
95  {
96  return i == 0 ? 1 : _ranks[ i - 1 ] * _diff( i - 1 );
97  }
98  void _inc( const size_type &m, const size_type &em, array< size_type > &i ) const
99  {
100  if( m != em && ++ i[ m ] == _ranks[ m ] )
101  {
102  i[ m ] = 0;
103  _inc( ( m + 1 ) % _mode, em, i );
104  }
105  }
106  void _out( ::std::ostream &o, const array< size_type > &d, array< size_type > &i, const size_type &m ) const
107  {
108  for( i[ m ] = 0 ; i[ m ] < _ranks[ m ] ; i[ m ] ++ )
109  {
110  if( m != 0 )
111  {
112  i[ m ] != 0 ? o << std::endl : o;
113  for( size_type ii = 0 ; ii < _mode - m - 1 ; o << " ", ii ++ );
114  m > 1 ? o << i[ m ] + 1 << " ----------" << ::std::endl : o;
115  _out( o, d, i, m - 1 );
116  }
117  else
118  {
119  o << operator [ ]( ::std::inner_product( &i[ 0 ], &i[ 0 ] + _mode, &d[ 0 ], 0 ) );
120  i[ m ] != _ranks[ m ] - 1 ? o << ", " : o;
121  }
122  }
123  }
124 
125  protected:
126  const size_type mode( ) const
127  {
128  return _mode;
129  }
130  const size_type rank( const size_type &m ) const
131  {
132  return _ranks[ m - 1 ];
133  }
134  const size_type size( ) const
135  {
136  return _data.size( );
137  }
138  value_type &operator [ ]( const size_type &i )
139  {
140  return _data[ i ];
141  }
142  const value_type &operator [ ]( const size_type &i ) const
143  {
144  return _data[ i ];
145  }
146  matrix_type unfold( const size_type &m ) const
147  {
148  array< size_type > i( _mode ), d( _mode );
149  for( size_type ii = 0 ; ii < _mode ; d[ ii ] = _diff( ii ), ii ++ );
150  matrix_type mat( _ranks[ m - 1 ], _data.size( ) / _ranks[ m - 1 ] );
151  for( typename matrix_type::size_type c = 0 ; c < mat.cols( ) ; c ++ )
152  {
153  for( typename matrix_type::size_type r = i[ m - 1 ] = 0 ; r < mat.rows( ) ; r ++, i[ m - 1 ] ++ )
154  {
155  mat( r, c ) = operator [ ]( ::std::inner_product( &i[ 0 ], &i[ 0 ] + i.size( ), &d[ 0 ], 0 ) );
156  }
157  _inc( m % _mode, m - 1, i );
158  }
159  return mat;
160  }
161  void out( ::std::ostream &o ) const
162  {
163  array< size_type > i( _mode ), d( _mode );
164  for( size_type ii = 0 ; ii < _mode ; d[ ii ] = _diff( ii ), ii ++ );
165  _out( o, d, i, _mode - 1 );
166  }
167 
168  _base_( ) : _mode( M ), _ranks( array< size_type >( M ) )
169  {
170  }
171  };
172 
173  template< int M, typename V, typename A >
174  class _data_ : public _base_< M, V, A >
175  {
176  _data_( )
177  {
178  }
179  };
180 
181  /*******************************************/
182  /* */
183  /* for 2nd order tensor */
184  /* */
185  /*******************************************/
186  template< typename V, typename A >
187  class _data_< 2, V, A >
188  {
189  template< int FM, typename FV, typename FA >
190  friend class tensor;
191 
192  typedef V value_type;
193  typedef A allocator;
194  typedef matrix< V, A > matrix_type;
195  typedef typename matrix_type::size_type size_type;
196 
197  matrix_type _mat;
198 
199  const size_type mode( ) const
200  {
201  return 2;
202  }
203  const size_type rank( const size_type &m ) const
204  {
205  return m == 1 ? _mat.rows( ) : _mat.cols( );
206  }
207  const size_type size( ) const
208  {
209  return _mat.size( );
210  }
211  void _resize( const size_type &r1, const size_type &r2 )
212  {
213  _mat.resize( r1, r2 );
214  }
215  void resize( const size_type &r1, const size_type &r2 )
216  {
217  _resize( r1, r2 );
218  _mat.fill( );
219  }
220  value_type &operator [ ]( const size_type &i )
221  {
222  return _mat[ i ];
223  }
224  const value_type &operator [ ]( const size_type &i ) const
225  {
226  return _mat[ i ];
227  }
228  value_type &operator ( )( const size_type &i1, const size_type &i2 )
229  {
230  return _mat( i1, i2 );
231  }
232  const value_type &operator ( )( const size_type &i1, const size_type &i2 ) const
233  {
234  return _mat( i1, i2 );
235  }
236  template< typename MV, typename MA >
237  _data_< 2, V, A > fold( const size_type &m, const matrix< MV, MA > &mat ) const
238  {
239  return _data_< 2, V, A >( m == 1 ? mat : mat.t( ) );
240  }
241  matrix_type unfold( const size_type &m ) const
242  {
243  return m == 1 ? _mat : _mat.t( );
244  }
245  void out( ::std::ostream &o ) const
246  {
247  o << _mat.t( );
248  }
249  template< typename AV, typename AA >
250  _data_< 2, V, A >( const size_t &r1, const size_t &r2, const array< AV, AA > &data ) : _mat( matrix_type( r1, r2 ) )
251  {
252  ::std::copy( &data[ 0 ], &data[ 0 ] + data.size( ), &_mat[ 0 ] );
253  }
254  template< typename AV, typename AA >
255  _data_< 2, V, A >( const array2< AV, AA > &img ) : _mat( matrix_type( img.width( ), img.height( ) ) )
256  {
257  ::std::copy( &img[ 0 ], &img[ 0 ] + img.size( ), &_mat[ 0 ] );
258  }
259  template< typename MV, typename MA >
260  _data_< 2, V, A >( const matrix< MV, MA > &mat ) : _mat( mat )
261  {
262  }
263  _data_< 2, V, A >( const size_type &r1, const size_type &r2 ) : _mat( matrix_type( r1, r2 ) )
264  {
265  }
266  _data_< 2, V, A >( ) : _mat( matrix_type( ) )
267  {
268  }
269  };
270 
271  /*******************************************/
272  /* */
273  /* for 3rd order tensor */
274  /* */
275  /*******************************************/
276  template< typename V, typename A >
277  class _data_< 3, V, A > : public _base_< 3, V, A >
278  {
279  template< int FM, typename FV, typename FA >
280  friend class tensor;
281 
282  typedef _base_< 3, V, A > base;
283  typedef typename base::value_type value_type;
284  typedef typename base::size_type size_type;
285 
286  void _resize( const size_type &r1, const size_type &r2, const size_type &r3 )
287  {
288  if( base::_ranks[ 0 ] != r1 || base::_ranks[ 1 ] != r2 || base::_ranks[ 2 ] != r3 )
289  {
290  base::_ranks[ 0 ] = r1;
291  base::_ranks[ 1 ] = r2;
292  base::_ranks[ 2 ] = r3;
293  base::_data.resize( r1 * r2 * r3 );
294  }
295  }
296  void resize( const size_type &r1, const size_type &r2, const size_type &r3 )
297  {
298  _resize( r1, r2, r3 );
299  base::_data.fill( );
300  }
301  value_type &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3 )
302  {
303  return base::_data[ i1 + i2 * base::_ranks[ 0 ] + i3 * base::_ranks[ 0 ] * base::_ranks[ 1 ] ];
304  }
305  const value_type &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3 ) const
306  {
307  return base::_data[ i1 + i2 * base::_ranks[ 0 ] + i3 * base::_ranks[ 0 ] * base::_ranks[ 1 ] ];
308  }
309  template< typename MV, typename MA >
310  _data_< 3, V, A > fold( const size_type &m, const matrix< MV, MA > &mat ) const
311  {
312  _data_< 3, V, A > ret( base::_ranks[ 0 ], base::_ranks[ 1 ], base::_ranks[ 2 ] );
313  ret._ranks[ m - 1 ] = mat.rows( );
314  ret._resize( base::_ranks[ 0 ], base::_ranks[ 1 ], base::_ranks[ 2 ] );
315  ret._fold( m, mat );
316  return ret;
317  }
318  template< typename AV, typename AA >
319  _data_< 3, V, A >( const size_type &r1, const size_type &r2, const size_type &r3, const array< AV, AA > &data )
320  {
321  _resize( r1, r2, r3 );
322  ::std::copy( &data[ 0 ], &data[ 0 ] + data.size( ), &base::_data[ 0 ] );
323  }
324  _data_< 3, V, A >( const size_type &r1, const size_type &r2, const size_type &r3 )
325  {
326  _resize( r1, r2, r3 );
327  }
328  template< typename TV, typename TA, typename AA >
329  _data_< 3, V, A >( const array< tensor< 2, TV, TA >, AA > &ts )
330  {
331  _resize( ts[ 0 ].rank( 1 ), ts[ 0 ].rank( 2 ), ts.size( ) );
332  for( size_type i = 0 ; i < ts.size( ) ; i ++ )
333  {
334  ::std::copy( &ts[ i ][ 0 ], &ts[ i ][ 0 ] + ts[ 0 ].size( ), &base::_data[ i * ts[ 0 ].size( ) ] );
335  }
336  }
337  _data_< 3, V, A >( )
338  {
339  }
340  };
341 
342  /*******************************************/
343  /* */
344  /* for 4th order tensor */
345  /* */
346  /*******************************************/
347  template< typename V, typename A >
348  class _data_< 4, V, A > : public _base_< 4, V, A >
349  {
350  template< int FM, typename FV, typename FA >
351  friend class tensor;
352 
353  typedef _base_< 4, V, A > base;
354  typedef typename base::value_type value_type;
355  typedef typename base::size_type size_type;
356 
357  void _resize( const size_type &r1, const size_type &r2, const size_type &r3, const size_type &r4 )
358  {
359  if( base::_ranks[ 0 ] != r1 || base::_ranks[ 1 ] != r2 || base::_ranks[ 2 ] != r3 || base::_ranks[ 3 ] != r4 )
360  {
361  base::_ranks[ 0 ] = r1;
362  base::_ranks[ 1 ] = r2;
363  base::_ranks[ 2 ] = r3;
364  base::_ranks[ 3 ] = r4;
365  base::_data.resize( r1 * r2 * r3 * r4 );
366  }
367  }
368  void resize( const size_type &r1, const size_type &r2, const size_type &r3, const size_type &r4 )
369  {
370  _resize( r1, r2, r3, r4 );
371  base::_data.fill( );
372  }
373  value_type &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3, const size_type &i4 )
374  {
375  return base::_data[ i1 + i2 * base::_ranks[ 0 ] + i3 * base::_ranks[ 0 ] * base::_ranks[ 1 ] + i4 * base::_ranks[ 0 ] * base::_ranks[ 1 ] * base::_ranks[ 2 ] ];
376  }
377  const value_type &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3, const size_type &i4 ) const
378  {
379  return base::_data[ i1 + i2 * base::_ranks[ 0 ] + i3 * base::_ranks[ 0 ] * base::_ranks[ 1 ] + i4 * base::_ranks[ 0 ] * base::_ranks[ 1 ] * base::_ranks[ 2 ] ];
380  }
381  template< typename MV, typename MA >
382  _data_< 4, V, A > fold( const size_type &m, const matrix< MV, MA > &mat ) const
383  {
384  _data_< 4, V, A > ret( base::_ranks[ 0 ], base::_ranks[ 1 ], base::_ranks[ 2 ], base::_ranks[ 3 ] );
385  ret.base::_ranks[ m - 1 ] = mat.rows( );
386  ret._resize( base::_ranks[ 0 ], base::_ranks[ 1 ], base::_ranks[ 2 ], base::_ranks[ 3 ] );
387  ret._fold( m, mat );
388  return ret;
389  }
390  template< typename AV, typename AA >
391  _data_< 4, V, A >( const size_type &r1, const size_type &r2, const size_type &r3, const size_type &r4, const array< AV, AA > &data )
392  {
393  _resize( r1, r2, r3, r4 );
394  ::std::copy( &data[ 0 ], &data[ 0 ] + data.size( ), &base::_data[ 0 ] );
395  }
396  _data_< 4, V, A >( const size_type &r1, const size_type &r2, const size_type &r3, const size_type &r4 )
397  {
398  _resize( r1, r2, r3, r4 );
399  }
400  template< typename TV, typename TA, typename AA >
401  _data_< 4, V, A >( const array< tensor< 3, TV, TA >, AA > &ts )
402  {
403  _resize( ts[ 0 ].rank( 1 ), ts[ 0 ].rank( 2 ), ts[ 0 ].rank( 3 ), ts.size( ) );
404  for( size_type i = 0 ; i < ts.size( ) ; i ++ )
405  {
406  ::std::copy( &ts[ i ][ 0 ], &ts[ i ][ 0 ] + ts[ 0 ].size( ), &base::_data[ i * ts[ 0 ].size( ) ] );
407  }
408  }
409  _data_< 4, V, A >( )
410  {
411  }
412  };
413 
414  /*******************************************/
415  /* */
416  /* for 5th order tensor */
417  /* */
418  /*******************************************/
419  template< typename V, typename A >
420  class _data_< 5, V, A > : public _base_< 5, V, A >
421  {
422  template< int FM, typename FV, typename FA >
423  friend class tensor;
424 
425  typedef _base_< 5, V, A > base;
426  typedef typename base::value_type value_type;
427  typedef typename base::size_type size_type;
428 
429  void _resize( const size_type &r1, const size_type &r2, const size_type &r3, const size_type &r4, const size_type &r5 )
430  {
431  if( base::_ranks[ 0 ] != r1 || base::_ranks[ 1 ] != r2 || base::_ranks[ 2 ] != r3 || base::_ranks[ 3 ] != r4 || base::_ranks[ 4 ] != r5 )
432  {
433  base::_ranks[ 0 ] = r1;
434  base::_ranks[ 1 ] = r2;
435  base::_ranks[ 2 ] = r3;
436  base::_ranks[ 3 ] = r4;
437  base::_ranks[ 4 ] = r5;
438  base::_data.resize( r1 * r2 * r3 * r4 * r5 );
439  }
440  }
441  void resize( const size_type &r1, const size_type &r2, const size_type &r3, const size_type &r4, const size_type &r5 )
442  {
443  _resize( r1, r2, r3, r4, r5 );
444  base::_data.fill( );
445  }
446  value_type &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3, const size_type &i4, const size_type &i5 )
447  {
448  return base::_data[ i1 + i2 * base::_ranks[ 0 ] + i3 * base::_ranks[ 0 ] * base::_ranks[ 1 ] + i4 * base::_ranks[ 0 ] * base::_ranks[ 1 ] * base::_ranks[ 2 ] + i5 * base::_ranks[ 0 ] * base::_ranks[ 1 ] * base::_ranks[ 2 ] * base::_ranks[ 3 ] ];
449  }
450  const value_type &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3, const size_type &i4, const size_type &i5 ) const
451  {
452  return base::_data[ i1 + i2 * base::_ranks[ 0 ] + i3 * base::_ranks[ 0 ] * base::_ranks[ 1 ] + i4 * base::_ranks[ 0 ] * base::_ranks[ 1 ] * base::_ranks[ 2 ] + i5 * base::_ranks[ 0 ] * base::_ranks[ 1 ] * base::_ranks[ 2 ] * base::_ranks[ 3 ] ];
453  }
454  template< typename MV, typename MA >
455  _data_< 5, V, A > fold( const size_type &m, const matrix< MV, MA > &mat ) const
456  {
457  _data_< 5, V, A > ret( base::_ranks[ 0 ], base::_ranks[ 1 ], base::_ranks[ 2 ], base::_ranks[ 3 ], base::_ranks[ 4 ] );
458  ret.base::_ranks[ m - 1 ] = mat.rows( );
459  ret._resize( base::_ranks[ 0 ], base::_ranks[ 1 ], base::_ranks[ 2 ], base::_ranks[ 3 ], base::_ranks[ 4 ] );
460  ret._fold( m, mat );
461  return ret;
462  }
463  template< typename AV, typename AA >
464  _data_< 5, V, A >( const size_type &r1, const size_type &r2, const size_type &r3, const size_type &r4, const size_type &r5, const array< AV, AA > &data )
465  {
466  _resize( r1, r2, r3, r4, r5 );
467  ::std::copy( &data[ 0 ], &data[ 0 ] + data.size( ), &base::_data[ 0 ] );
468  }
469  _data_< 5, V, A >( const size_type &r1, const size_type &r2, const size_type &r3, const size_type &r4, const size_type &r5 )
470  {
471  _resize( r1, r2, r3, r4, r5 );
472  }
473  template< typename TV, typename TA, typename AA >
474  _data_< 5, V, A >( const array< tensor< 4, TV, TA >, AA > &ts )
475  {
476  _resize( ts[ 0 ].rank( 1 ), ts[ 0 ].rank( 2 ), ts[ 0 ].rank( 3 ), ts[ 0 ].rank( 4 ), ts.size( ) );
477  for( size_type i = 0 ; i < ts.size( ) ; i ++ )
478  {
479  ::std::copy( &ts[ i ][ 0 ], &ts[ i ][ 0 ] + ts[ 0 ].size( ), &base::_data[ i * ts[ 0 ].size( ) ] );
480  }
481  }
482  _data_< 5, V, A >( )
483  {
484  }
485  };
486 }
487 
497 template< int M, typename V, typename A = ::std::allocator< V > >
498 class tensor
499 {
500  typedef matrix< V, A > matrix_type;
501  typedef typename __tensor__::_data_< M, V, A > data_type;
502 
503 public:
504  typedef V value_type;
505  typedef A allocator;
506  typedef typename data_type::size_type size_type;
507 
508 private:
509  data_type _d;
510 
511  bool _is_valid( const size_type &m ) const
512  {
513  return 0 < m && m <= mode( );
514  }
515  template< typename TV, typename TA >
516  bool _is_same_size( const tensor< M, TV, TA > &t ) const
517  {
518  for( size_type i = 1 ; i <= mode( ) ; i ++ )
519  {
520  if( rank( i ) != t.rank( i ) )
521  {
522  return false;
523  }
524  }
525  return true;
526  }
527  template< typename MV, typename MA >
528  tensor< M, V, A > _fold( const size_type &m, const matrix< MV, MA > &mat ) const
529  {
530  tensor< M, V, A > ret;
531  ret._d = _d.fold( m, mat );
532  return ret;
533  }
534  matrix_type _unfold( const size_type &m ) const
535  {
536  return _d.unfold( m );
537  }
538  template< typename MV, typename MA >
539  tensor< M, V, A > _x( const size_type &m, const matrix< MV, MA > &mat ) const
540  {
541  return _fold( m, mat * _unfold( m ) );
542  }
543 
544 public:
546  const size_type mode( ) const
547  {
548  return _d.mode( );
549  }
555  const size_type rank( const size_type &mode ) const
556  {
557  if( _is_valid( mode ) )
558  {
559  return _d.rank( mode );
560  }
561  else
562  {
563  ::std::cerr << "can't obtain rank" << ::std::endl;
564  return 0;
565  }
566  }
568  const size_type size( ) const
569  {
570  return _d.size( );
571  }
577  value_type &operator [ ]( const size_type &i )
578  {
579  return _d[ i ];
580  }
586  const value_type &operator [ ]( const size_type &i ) const
587  {
588  return _d[ i ];
589  }
600  matrix_type unfold( const size_type &mode ) const
601  {
602  if( _is_valid( mode ) )
603  {
604  return _unfold( mode );
605  }
606  else
607  {
608  ::std::cerr << "can't unfold tensor" << ::std::endl;
609  return matrix_type( );
610  }
611  }
620  template< typename MV, typename MA >
621  tensor< M, V, A > x( const size_type &mode, const matrix< MV, MA > &mat ) const
622  {
623  if( _is_valid( mode ) )
624  {
625  return _x( mode, mat );
626  }
627  else
628  {
629  ::std::cerr << "can't obtain n-mode product" << ::std::endl;
630  return tensor< M, V, A >( );
631  }
632  }
633  void out( ::std::ostream &o ) const
634  {
635  _d.out( o );
636  }
638  template< typename TV, typename TA >
639  bool operator ==( const tensor< M, TV, TA > &t ) const
640  {
641  if( !_is_same_size( t ) )
642  {
643  return false;
644  }
645  else
646  {
647  for( size_type i = 0 ; i < size( ) ; i ++ )
648  {
649  if( operator [ ]( i ) != t[ i ] )
650  {
651  return false;
652  }
653  }
654  return true;
655  }
656  }
658  template< typename TV, typename TA >
659  bool operator !=( const tensor< M, TV, TA > &t ) const
660  {
661  return !operator ==( t );
662  }
664  template< typename TV, typename TA >
665  const tensor< M, V, A > &operator +=( const tensor< M, TV, TA > &t )
666  {
667  if( _is_same_size( t ) )
668  {
669  ::std::transform( &operator [ ]( 0 ), &operator [ ]( 0 ) + size( ), &t[ 0 ], &operator [ ]( 0 ), ::std::plus< value_type >( ) );
670  }
671  else
672  {
673  ::std::cerr << "can't add tensors" << ::std::endl;
674  }
675  return *this;
676  }
678  template< typename TV, typename TA >
679  const tensor< M, V, A > &operator -=( const tensor< M, TV, TA > &t )
680  {
681  if( _is_same_size( t ) )
682  {
683  ::std::transform( &operator [ ]( 0 ), &operator [ ]( 0 ) + size( ), &t[ 0 ], &operator [ ]( 0 ), ::std::minus< value_type >( ) );
684  }
685  else
686  {
687  ::std::cerr << "can't subtract tensors" << ::std::endl;
688  }
689  return *this;
690  }
692  template< typename VV >
693  const tensor< M, V, A > &operator *=( const VV &v )
694  {
695  ::std::transform( &operator [ ]( 0 ), &operator [ ]( 0 ) + size( ), &operator [ ]( 0 ), ::std::bind2nd( ::std::multiplies< value_type >( ), v ) );
696  return *this;
697  }
699  template< typename VV >
700  const tensor< M, V, A > &operator /=( const VV &v )
701  {
702  ::std::transform( &operator [ ]( 0 ), &operator [ ]( 0 ) + size( ), &operator [ ]( 0 ), ::std::bind2nd( ::std::divides< value_type >( ), v ) );
703  return *this;
704  }
706  tensor< M, V, A >( ) : _d( )
707  {
708  }
709 
710 
711  /*******************************************/
712  /* */
713  /* for 2nd order tensor */
714  /* */
715  /*******************************************/
716 
722  void resize( const size_type &rank1, const size_type &rank2 )
723  {
724  _d.resize( rank1, rank2 );
725  }
732  value_type &operator ( )( const size_type &i1, const size_type &i2 )
733  {
734  return _d( i1, i2 );
735  }
742  const value_type &operator ( )( const size_type &i1, const size_type &i2 ) const
743  {
744  return _d( i1, i2 );
745  }
757  template< typename MV, typename MA >
758  void hosvd( tensor< M, V, A > &z, matrix< MV, MA > &u1, matrix< MV, MA > &u2 ) const
759  {
760  matrix< MV, MA > s;
761  svd( _unfold( 1 ), u1, s, u2 );
762  z = _x( 1, u1.dagger( ) )._x( 2, u2 );
763  u2 = u2.dagger( );
764  }
771  template< typename AV, typename AA >
772  tensor< M, V, A >( const size_type &rank1, const size_type &rank2, const array< AV, AA > &data ) : _d( rank1, rank2, data )
773  {
774  }
781  template< typename AV, typename AA >
782  tensor< M, V, A >( const array2< AV, AA > &img ) : _d( img.width( ), img.height( ), img )
783  {
784  }
791  template< typename MV, typename MA >
792  tensor< M, V, A >( const matrix< MV, MA > &mat ) : _d( mat )
793  {
794  }
800  tensor< M, V, A >( const size_type &rank1, const size_type &rank2 ) : _d( rank1, rank2 )
801  {
802  }
803 
804 
805  /*******************************************/
806  /* */
807  /* for 3rd order tensor */
808  /* */
809  /*******************************************/
810 
817  void resize( const size_type &rank1, const size_type &rank2, const size_type &rank3 )
818  {
819  _d.resize( rank1, rank2, rank3 );
820  }
828  value_type &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3 )
829  {
830  return _d( i1, i2, i3 );
831  }
839  const value_type &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3 ) const
840  {
841  return _d( i1, i2, i3 );
842  }
853  template< typename MV, typename MA >
854  void hosvd( tensor< M, V, A > &z, matrix< MV, MA > &u1, matrix< MV, MA > &u2, matrix< MV, MA > &u3 ) const
855  {
856  matrix< MV, MA > s, tmp;
857  svd( _unfold( 1 ), u1, s, tmp );
858  svd( _unfold( 2 ), u2, s, tmp );
859  svd( _unfold( 3 ), u3, s, tmp );
860  z = _x( 1, u1.dagger( ) )._x( 2, u2.dagger( ) )._x( 3, u3.dagger( ) );
861  }
869  template< typename AV, typename AA >
870  tensor< M, V, A >( const size_type &rank1, const size_type &rank2, const size_type &rank3, const array< AV, AA > &data ) : _d( rank1, rank2, rank3, data )
871  {
872  }
879  tensor< M, V, A >( const size_type &rank1, const size_type &rank2, const size_type &rank3 ) : _d( rank1, rank2, rank3 )
880  {
881  }
886  template< typename TV, typename TA, typename AA >
887  tensor< M, V, A >( const array< tensor< 2, TV, TA >, AA > &ts ) : _d( ts )
888  {
889  }
890 
891 
892  /*******************************************/
893  /* */
894  /* for 4th order tensor */
895  /* */
896  /*******************************************/
897 
905  void resize( const size_type &rank1, const size_type &rank2, const size_type &rank3, const size_type &rank4 )
906  {
907  _d.resize( rank1, rank2, rank3, rank4 );
908  }
917  value_type &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3, const size_type &i4 )
918  {
919  return _d( i1, i2, i3, i4 );
920  }
929  const value_type &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3, const size_type &i4 ) const
930  {
931  return _d( i1, i2, i3, i4 );
932  }
944  template< typename MV, typename MA >
945  void hosvd( tensor< M, V, A > &z, matrix< MV, MA > &u1, matrix< MV, MA > &u2, matrix< MV, MA > &u3, matrix< MV, MA > &u4 ) const
946  {
947  matrix< MV, MA > s, tmp;
948  svd( _unfold( 1 ), u1, s, tmp );
949  svd( _unfold( 2 ), u2, s, tmp );
950  svd( _unfold( 3 ), u3, s, tmp );
951  svd( _unfold( 4 ), u4, s, tmp );
952  z = _x( 1, u1.dagger( ) )._x( 2, u2.dagger( ) )._x( 3, u3.dagger( ) )._x( 4, u4.dagger( ) );
953  }
962  template< typename AV, typename AA >
963  tensor< M, V, A >( const size_type &rank1, const size_type &rank2, const size_type &rank3, const size_type &rank4, const array< AV, AA > &data ) : _d( rank1, rank2, rank3, rank4, data )
964  {
965  }
973  tensor< M, V, A >( const size_type &rank1, const size_type &rank2, const size_type &rank3, const size_type &rank4 ) : _d( rank1, rank2, rank3, rank4 )
974  {
975  }
980  template< typename TV, typename TA, typename AA >
981  tensor< M, V, A >( const array< tensor< 3, TV, TA >, AA > &ts ) : _d( ts )
982  {
983  }
984 
985 
986  /*******************************************/
987  /* */
988  /* for 5th order tensor */
989  /* */
990  /*******************************************/
991 
1000  void resize( const size_type &rank1, const size_type &rank2, const size_type &rank3, const size_type &rank4, const size_type &rank5 )
1001  {
1002  _d.resize( rank1, rank2, rank3, rank4, rank5 );
1003  }
1013  tensor< M, V, A > &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3, const size_type &i4, const size_type &i5 )
1014  {
1015  return _d( i1, i2, i3, i4, i5 );
1016  }
1026  const tensor< M, V, A > &operator ( )( const size_type &i1, const size_type &i2, const size_type &i3, const size_type &i4, const size_type &i5 ) const
1027  {
1028  return _d( i1, i2, i3, i4, i5 );
1029  }
1042  template< typename MV, typename MA >
1043  void hosvd( tensor< M, V, A > &z, matrix< MV, MA > &u1, matrix< MV, MA > &u2, matrix< MV, MA > &u3, matrix< MV, MA > &u4, matrix< MV, MA > &u5 ) const
1044  {
1045  matrix< MV, MA > s, tmp;
1046  svd( _unfold( 1 ), u1, s, tmp );
1047  svd( _unfold( 2 ), u2, s, tmp );
1048  svd( _unfold( 3 ), u3, s, tmp );
1049  svd( _unfold( 4 ), u4, s, tmp );
1050  svd( _unfold( 5 ), u5, s, tmp );
1051  z = _x( 1, u1.dagger( ) )._x( 2, u2.dagger( ) )._x( 3, u3.dagger( ) )._x( 4, u4.dagger( ) )._x( 5, u5.dagger( ) );
1052  }
1062  template< typename AV, typename AA >
1063  tensor< M, V, A >( const size_type &rank1, const size_type &rank2, const size_type &rank3, const size_type &rank4, const size_type &rank5, const array< AV, AA > &data ) : _d( rank1, rank2, rank3, rank4, rank5, data )
1064  {
1065  }
1074  tensor< M, V, A >( const size_type &rank1, const size_type &rank2, const size_type &rank3, const size_type &rank4, const size_type &rank5 ) : _d( rank1, rank2, rank3, rank4, rank5 )
1075  {
1076  }
1081  template< typename TV, typename TA, typename AA >
1082  tensor< M, V, A >( const array< tensor< 4, TV, TA >, AA > &ts ) : _d( ts )
1083  {
1084  }
1085 };
1086 
1087 
1094 template< int M, typename V, typename A >
1095 inline tensor< M, V, A > operator +( const tensor< M, V, A > &t1, const tensor< M, V, A > &t2 )
1096 {
1097  return tensor< M, V, A >( t1 ) += t2;
1098 }
1099 
1106 template< int M, typename V, typename A >
1107 inline tensor< M, V, A > operator -( const tensor< M, V, A > &t1, const tensor< M, V, A > &t2 )
1108 {
1109  return tensor< M, V, A >( t1 ) -= t2;
1110 }
1111 
1118 template< int M, typename V, typename A >
1119 inline tensor< M, V, A > operator *( const tensor< M, V, A > &t, const typename type_trait< V >::value_type &v )
1120 {
1121  return tensor< M, V, A >( t ) *= v;
1122 }
1123 
1130 template< int M, typename V, typename A >
1131 inline tensor< M, V, A > operator *( const typename type_trait< V >::value_type &v, const tensor< M, V, A > &t )
1132 {
1133  return tensor< M, V, A >( t ) *= v;
1134 }
1135 
1142 template< int M, typename V, typename A >
1143 inline tensor< M, V, A > operator /( const tensor< M, V, A > &t, const typename type_trait< V >::value_type &v )
1144 {
1145  return tensor< M, V, A >( t ) /= v;
1146 }
1147 
1154 template < int M, typename V, typename A >
1155 inline ::std::ostream &operator <<( ::std::ostream &o, const tensor< M, V, A > &t )
1156 {
1157  t.out( o );
1158  return o;
1159 }
1160 
1161 
1162  /*******************************************/
1163  /* */
1164  /* for 2nd order tensor */
1165  /* */
1166  /*******************************************/
1167 
1180 template< typename TV, typename TA, typename MV, typename MA >
1181 inline void hosvd( const tensor< 2, TV, TA > &t, tensor< 2, TV, TA > &z, matrix< MV, MA > &u1, matrix< MV, MA > &u2 )
1182 {
1183  t.hosvd( z, u1, u2 );
1184 }
1185 
1186 
1187  /*******************************************/
1188  /* */
1189  /* for 3rd order tensor */
1190  /* */
1191  /*******************************************/
1192 
1204 template< typename TV, typename TA, typename MV, typename MA >
1205 inline void hosvd( const tensor< 3, TV, TA > &t, tensor< 3, TV, TA > &z, matrix< MV, MA > &u1, matrix< MV, MA > &u2, matrix< MV, MA > &u3 )
1206 {
1207  t.hosvd( z, u1, u2, u3 );
1208 }
1209 
1210 
1211  /*******************************************/
1212  /* */
1213  /* for 4th order tensor */
1214  /* */
1215  /*******************************************/
1216 
1229 template< typename TV, typename TA, typename MV, typename MA >
1230 inline void hosvd( const tensor< 4, TV, TA > &t, tensor< 4, TV, TA > &z, matrix< MV, MA > &u1, matrix< MV, MA > &u2, matrix< MV, MA > &u3, matrix< MV, MA > &u4 )
1231 {
1232  t.hosvd( z, u1, u2, u3, u4 );
1233 }
1234 
1235 
1236  /*******************************************/
1237  /* */
1238  /* for 5th order tensor */
1239  /* */
1240  /*******************************************/
1241 
1255 template< typename TV, typename TA, typename MV, typename MA >
1256 inline void hosvd( const tensor< 5, TV, TA > &t, tensor< 5, TV, TA > &z, matrix< MV, MA > &u1, matrix< MV, MA > &u2, matrix< MV, MA > &u3, matrix< MV, MA > &u4, matrix< MV, MA > &u5 )
1257 {
1258  t.hosvd( z, u1, u2, u3, u4, u5 );
1259 }
1260 
1261 
1262 // closing of mist namespace
1263 _MIST_END
1264 
1265 #endif // __INCLUDE_MIST_TENSOR__
Generated on Wed Nov 12 2014 19:44:25 for MIST by doxygen 1.8.1.2