- メインページ
- モジュール
- ネームスペース
- データ構造
- ファイル

- ファイル一覧
- グローバル

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__