vmmlib  1.7.0
 All Classes Namespaces Functions Pages
tensor4.hpp
1 /*
2  * Copyright (c) 2006-2014, Visualization and Multimedia Lab,
3  * University of Zurich <http://vmml.ifi.uzh.ch>,
4  * Eyescale Software GmbH,
5  * Blue Brain Project, EPFL
6  *
7  * This file is part of VMMLib <https://github.com/VMML/vmmlib/>
8  *
9  * Redistribution and use in source and binary forms, with or without
10  * modification, are permitted provided that the following conditions are met:
11  *
12  * Redistributions of source code must retain the above copyright notice, this
13  * list of conditions and the following disclaimer. Redistributions in binary
14  * form must reproduce the above copyright notice, this list of conditions and
15  * the following disclaimer in the documentation and/or other materials provided
16  * with the distribution. Neither the name of the Visualization and Multimedia
17  * Lab, University of Zurich nor the names of its contributors may be used to
18  * endorse or promote products derived from this software without specific prior
19  * written permission.
20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
21  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
24  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
25  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
26  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
27  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
28  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
29  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30  * POSSIBILITY OF SUCH DAMAGE.
31  */
32 
33 /* @author David Klaper
34  * @author Susanne Suter
35  * @author Jonas Boesch
36  *
37  * a tensor is a generalization of a multidimensional array
38  * a tensor4 is a tensor data structure with four modes I1, I2, I3 and I4
39  */
40 
41 #ifndef __VMML__TENSOR4__HPP__
42 #define __VMML__TENSOR4__HPP__
43 
44 #include <fstream> // file I/O
45 #include <vmmlib/enable_if.hpp>
46 #include "tensor3.hpp"
47 
48 
49 namespace vmml
50 {
51 
52  // tensor with four modes, containing a series I4 x tensor4 as (I3 of I1 x
53  // I2 vmml matrices). I1 is number of rows, I2 is number of columns and I3
54  // is number of tubes, I4 is number of tensor4s
55 
56  template< size_t I1, size_t I2, size_t I3, size_t I4, typename T = float >
57  class tensor4
58  {
59  public:
60  typedef T value_type;
61  typedef T* pointer;
62  typedef T& reference;
63 
65 
66  //TODO: maybe tensor4 iterator
67  //TODO: unfolding along all modes
68  //TODO: accessors to tensor3 (along all modes)
69 
70  // Average all values along 1 axis
71  tensor3_t& average_I4(tensor3_t& t3) const;
72 
77 
78  inline void get_tensor3( const size_t i4_, tensor3_t& t3_data_ ) const;
79 
80  inline tensor3_t& get_tensor3( size_t i4_ );
81  inline const tensor3_t& get_tensor3( size_t i4_ ) const;
82  inline void set_tensor3( size_t i4_, const tensor3_t& t3_data_ );
83 
84  static const size_t ROWS = I1;
85  static const size_t COLS = I2;
86  static const size_t SLICES = I3;
87  static const size_t T3S = I4;
88  static const size_t MATRIX_SIZE = I1 * I2;
89  static const size_t T3_SIZE = I1 * I2 * I3;
90  static const size_t SIZE = I1 * I2 * I3 * I4;
91 
92  static size_t get_array_size_in_bytes();
93 
94  // WARNING: dangerous. Use before destruction if you want to prevent
95  // a delete call for the assigned T* _array in the destructor.
96  void clear_array_pointer();
97 
98  // accessors
99  inline T& operator()( size_t i1, size_t i2, size_t i3, size_t i4 );
100  inline const T& operator()( size_t i1, size_t i2, size_t i3, size_t i4 ) const;
101 
102  inline T& at( size_t i1, size_t i2, size_t i3, size_t i4 );
103  inline const T& at( size_t i1, size_t i2, size_t i3, size_t i4 ) const;
104 
105 
106  // ctors
107  tensor4();
108 
109  explicit tensor4(void* memory) ;
110 
111  tensor4( const tensor4& source );
112 
113  template< typename U >
114  tensor4( const tensor4< I1, I2, I3, I4, U >& source_ );
115 
116  template< size_t J1, size_t J2, size_t J3, size_t J4 >
117  tensor4( const tensor4< J1, J2, J3, J4, T >& source_ );
118 
119  ~tensor4();
120 
121  size_t size() const; // return I1 * I2 * I3 * I4;
122 
123  // sets all elements to fill_value
124  void operator=( T fill_value );
125  void fill( T fill_value ); //special case of set method (all values are set to the same value!)
126 
127  //sets all tensor values with random values
128  //if seed is negative, srand( seed ) should have been set outside fill_random
129  //if seed is 0 or greater srand( seed ) will be called with the given seed
130  void fill_random( int seed = -1 );
131  void fill_random_signed( int seed = -1 );
132  void fill_increasing_values( );
133 
134  const tensor4& operator=( const tensor4& source_ );
135 
136 
137  // note: this function copies elements until either the matrix is full or
138  // the iterator equals end_.
139  template< typename input_iterator_t >
140  void set( input_iterator_t begin_, input_iterator_t end_,
141  bool row_major_layout = true );
142  void zero();
143 
144  T get_min() const;
145  T get_max() const;
146  T get_abs_min() const;
147  T get_abs_max() const;
148 
149  //returns number of non-zeros
150  size_t nnz() const;
151  size_t nnz( const T& threshold_ ) const;
152  void threshold( const T& threshold_value_ );
153 
154  //error computation
155  double frobenius_norm() const;
156  double frobenius_norm( const tensor4< I1, I2, I3, I4, T >& other ) const;
157  double avg_frobenius_norm() const;
158  double mse( const tensor4< I1, I2, I3, I4, T >& other ) const; // mean-squared error
159  double rmse( const tensor4< I1, I2, I3, I4, T >& other ) const; //root mean-squared error
160  double compute_psnr( const tensor4< I1, I2, I3, I4, T >& other, const T& max_value_ ) const; //peak signal-to-noise ratio
161  void mean( T& mean_ ) const;
162  double mean() const;
163  double variance() const;
164  double stdev() const;
165 
166  template< typename TT >
167  void cast_from( const tensor4< I1, I2, I3, I4, TT >& other );
168 
169  template< size_t J1, size_t J2, size_t J3, size_t J4, typename TT >
170  void cast_from( const tensor4< J1, J2, J3, J4, TT >& other );
171 
172 
173  template< typename TT >
174  void float_t_to_uint_t( const tensor4< I1, I2, I3, I4, TT >& other );
175 
176  //check if corresponding tensor values are equal or not
177  bool operator==( const tensor4& other ) const;
178  bool operator!=( const tensor4& other ) const;
179 
180  // due to limited precision, two 'idential' tensor4 might seem different.
181  // this function allows to specify a tolerance when comparing matrices.
182  bool equals( const tensor4& other, T tolerance ) const;
183  // this version takes a comparison functor to compare the components of
184  // the two tensor4 data structures
185  template< typename compare_t >
186  bool equals( const tensor4& other, compare_t& cmp ) const;
187 
188 
189  inline tensor4 operator+( T scalar ) const;
190  inline tensor4 operator-( T scalar ) const;
191 
192  void operator+=( T scalar );
193  void operator-=( T scalar );
194 
195  inline tensor4 operator+( const tensor4& other ) const;
196  inline tensor4 operator-( const tensor4& other ) const;
197 
198  template< size_t J1, size_t J2, size_t J3, size_t J4 >
200  operator+=( const tensor4< J1, J2, J3, J4, T>& other );
201 
202  void operator+=( const tensor4& other );
203  void operator-=( const tensor4& other );
204 
205  //
206  // tensor4-scalar operations / scaling
207  //
208  tensor4 operator*( T scalar );
209  void operator*=( T scalar );
210 
211  tensor4 operator/( T scalar );
212  void operator/=( T scalar );
213 
214 
215  inline tensor4< I1, I2, I3, I4, T > operator-() const;
216  tensor4< I1, I2, I3, I4, T > negate() const;
217 
218  friend std::ostream& operator << ( std::ostream& os,
219  const tensor4< I1,I2,I3,I4,T >& t4 )
220  {
221  //FIXME: to this directly with tensors
222  //sth like:
223  //os << t4.get_tensor3( i ) << " xxx " << std::endl;
224  for(size_t i4 = 0; i4 < I4; ++i4)
225  {
226  for(size_t i3 = 0; i3 < I3; ++i3)
227  {
228  for(size_t i1 = 0; i1 < I1; ++i1)
229  {
230  os << "(";
231 
232  for(size_t i2 = 0; i2 < I2; ++i2)
233  {
234  if( i2+1 == I2 )
235  {
236  os << T(t4.at( i1, i2, i3, i4 )) ;
237  } else {
238  os << T(t4.at( i1, i2, i3, i4 )) << ", " ;
239  }
240 
241  }
242  os << ")" << std::endl;
243  }
244  if ( i3+1 < I3 )
245  os << " *** " << std::endl;
246  }
247  if ( i4 + 1 < I4 )
248  {
249  os << "---- " << std::endl;
250  }
251  }
252  return os;
253  }
254 
255  template< size_t J1, size_t J2, size_t J3, size_t J4 >
257  get_sub_tensor4(tensor4<J1, J2, J3, J4, T >& result, size_t row_offset, size_t col_offset, size_t slice_offset, size_t t3_offset) const;
258 
259  // static members
260  static void tensor4_allocate_data( T*& array_ );
261  static void tensor4_deallocate_data( T*& array_ );
262 
263  static const tensor4< I1, I2, I3, I4, T > ZERO;
264 
265  T* get_array_ptr();
266  const T* get_array_ptr() const;
267 
268  void mode1_unfolding_fwd(mode1_unfolding_type& unfolding) const;
269  void mode2_unfolding_fwd(mode2_unfolding_type& unfolding) const;
270  void mode3_unfolding_fwd(mode3_unfolding_type& unfolding) const;
271  void mode4_unfolding_fwd(mode4_unfolding_type& unfolding) const;
272 
273  // computes the array index for direct access
274  inline size_t compute_index( size_t i1, size_t i2, size_t i3, size_t i4 ) const;
275 
276  template< typename TT >
277  void quantize_log(tensor4< I1, I2, I3, I4, TT >& quantized_, tensor4< I1, I2, I3, I4, char >& signs_, T& min_value_, T& max_value_, const TT& tt_range_) const;
278 
279  protected:
280  tensor3_t& _get_tensor3( size_t index_ );
281  const tensor3_t& _get_tensor3( size_t index_ ) const;
282 
283  T* _array;
284 
285  }; // class tensor4
286 
287 #define VMML_TEMPLATE_STRING template< size_t I1, size_t I2, size_t I3, size_t I4, typename T >
288 #define VMML_TEMPLATE_CLASSNAME tensor4< I1, I2, I3, I4, T >
289 
290 
291 
292  // WARNING: make sure that memory is a pointer to a memory block of
293  // sufficient size (that is, is at least get_array_size_in_bytes())
294  VMML_TEMPLATE_STRING
295  VMML_TEMPLATE_CLASSNAME::tensor4( void* memory )
296  : _array( reinterpret_cast<T*>( memory ) )
297  {
298  assert( _array );
299  }
300 
301 
302  VMML_TEMPLATE_STRING
303  VMML_TEMPLATE_CLASSNAME::tensor4()
304  : _array()
305  {
306  tensor4_allocate_data( _array );
307  }
308 
309  VMML_TEMPLATE_STRING
310  VMML_TEMPLATE_CLASSNAME::~tensor4()
311  {
312  tensor4_deallocate_data( _array );
313  }
314 
315  VMML_TEMPLATE_STRING
316  template< size_t J1, size_t J2, size_t J3, size_t J4 >
317  typename enable_if< J1 <= I1 && J2 <= I2 && J3 <= I3 && J4 <= I4 >::type*
318  VMML_TEMPLATE_CLASSNAME::
319  get_sub_tensor4(tensor4<J1, J2, J3, J4, T >& result,
320  size_t row_offset, size_t col_offset, size_t slice_offset, size_t t3_offset) const {
321  #ifdef VMMLIB_SAFE_ACCESSORS
322  if (J1 + row_offset > I1 || J2 + col_offset > I2 || J3 + slice_offset > I3 || J4 + t3_offset > I4)
323  VMMLIB_ERROR("index out of bounds.", VMMLIB_HERE);
324  #endif
325 
326  for (size_t t3 = 0; t3 < J4; ++t3) {
327  for (size_t slice = 0; slice < J3; ++slice) {
328  for (size_t row = 0; row < J1; ++row) {
329  for (size_t col = 0; col < J2; ++col) {
330  result.at(row, col, slice, t3)
331  = at(row_offset + row, col_offset + col, slice_offset + slice, t3_offset + t3);
332  }
333  }
334  }
335  }
336  return 0;
337  }
338 
339  VMML_TEMPLATE_STRING
340  void
341  VMML_TEMPLATE_CLASSNAME::
342  tensor4_allocate_data( T*& array_ )
343  {
344  array_ = new T[ I1 * I2 * I3 * I4];
345  }
346 
347  VMML_TEMPLATE_STRING
348  void
349  VMML_TEMPLATE_CLASSNAME::
350  tensor4_deallocate_data( T*& array_ )
351  {
352  if (array_)
353  {
354  delete[] array_;
355  }
356  }
357 
358 
359 
360  VMML_TEMPLATE_STRING
361  T*
362  VMML_TEMPLATE_CLASSNAME::get_array_ptr()
363  {
364  return _array;
365  }
366 
367 
368 
369  VMML_TEMPLATE_STRING
370  const T*
371  VMML_TEMPLATE_CLASSNAME::get_array_ptr() const
372  {
373  return _array;
374  }
375 
376  VMML_TEMPLATE_STRING
377  VMML_TEMPLATE_CLASSNAME::tensor4( const tensor4& source_ )
378  {
379  tensor4_allocate_data( _array );
380  (*this) = source_;
381 
382  }
383 
384  VMML_TEMPLATE_STRING
385  template< typename U >
386  VMML_TEMPLATE_CLASSNAME::tensor4( const tensor4< I1, I2, I3, I4, U >& source_ )
387  {
388  const U* s_array = source_.get_array_ptr();
389  tensor4_allocate_data( _array );
390  for (size_t index = 0; index < SIZE; ++index)
391  {
392  _array[ index ] = static_cast< T >( s_array[ index ] );
393  }
394  }
395 
396  VMML_TEMPLATE_STRING
397  template< size_t J1, size_t J2, size_t J3, size_t J4 >
398  VMML_TEMPLATE_CLASSNAME::tensor4( const tensor4< J1, J2, J3, J4, T >& source_ )
399  {
400  tensor4_allocate_data( _array );
401  size_t min1 = J1 < I1 ? J1 : I1;
402  size_t min2 = J2 < I2 ? J2 : I2;
403  size_t min3 = J3 < I3 ? J3 : I3;
404  size_t min4 = J4 < I4 ? J4 : I4;
405 
406  zero();
407 
408  for(size_t i4 = 0; i4 < min4; ++i4 )
409  {
410  for(size_t i3 = 0; i3 < min3; ++i3)
411  {
412  for(size_t i2 = 0; i2 < min2; ++i2)
413  {
414  for(size_t i1 = 0; i1 < min1; ++i1)
415  {
416  at(i1, i2, i3, i4) = source_(i1, i2, i3, i4);
417  }
418  }
419  }
420  }
421 
422  }
423 
424  VMML_TEMPLATE_STRING
425  template< typename input_iterator_t >
426  void
427  VMML_TEMPLATE_CLASSNAME::set(input_iterator_t begin_, input_iterator_t end_, bool row_major_layout) {
428  input_iterator_t it(begin_);
429  if (row_major_layout) {
430  for (size_t i4 = 0; i4 < I4; ++i4) {
431  for (size_t i3 = 0; i3 < I3; ++i3) {
432  for (size_t i1 = 0; i1 < I1; ++i1) {
433  for (size_t i2 = 0; i2 < I2; ++i2, ++it) {
434  if (it == end_)
435  return;
436  at(i1, i2, i3, i4) = static_cast<T> (*it);
437  }
438  }
439  }
440  }
441  } else {
442  VMMLIB_ERROR( "Tensor4: set() not implemented for non-row major", VMMLIB_HERE );
443  // TODO
444 // std::copy(it, it + (I1 * I2 * I3), begin());
445  }
446  }
447 
448  VMML_TEMPLATE_STRING
449  bool
450  VMML_TEMPLATE_CLASSNAME::equals(const tensor4< I1, I2, I3, I4, T >& other, T tolerance) const {
451  bool is_ok = true;
452  if (T3S != other.T3S) {
453  return false;
454  }
455  for (size_t i4 = 0; is_ok && i4 < I4; ++i4) {
456  is_ok = _get_tensor3(i4).equals(other._get_tensor3(i4), tolerance);
457  }
458  return is_ok;
459  }
460 
461  VMML_TEMPLATE_STRING
462  size_t
463  VMML_TEMPLATE_CLASSNAME::size() const
464  {
465  return SIZE;
466  }
467 
468  VMML_TEMPLATE_STRING
469  bool
470  VMML_TEMPLATE_CLASSNAME::operator==( const tensor4< I1, I2, I3, I4, T >& other ) const
471  {
472  const T* other_array = other.get_array_ptr();
473  for(size_t index = 0; index < SIZE; ++index)
474  {
475  if(_array[ index ] != other_array[ index] )
476  {
477  return false;
478  }
479  }
480 
481  return true;
482  }
483 
484 
485  VMML_TEMPLATE_STRING
486  bool
487  VMML_TEMPLATE_CLASSNAME::operator!=( const tensor4< I1, I2, I3, I4, T >& other ) const
488  {
489  return ! operator==( other );
490  }
491 
492 
493  VMML_TEMPLATE_STRING
494  typename VMML_TEMPLATE_CLASSNAME::tensor3_t&
495  VMML_TEMPLATE_CLASSNAME::
496  _get_tensor3( size_t index_ )
497  {
498  tensor3<I1, I2, I3, T>* tens3 = new tensor3<I1, I2, I3, T>( (void *)(_array + I1 * I2 * I3 * index_));
499  return *tens3;
500  }
501 
502 
503  VMML_TEMPLATE_STRING
504  const typename VMML_TEMPLATE_CLASSNAME::tensor3_t&
505  VMML_TEMPLATE_CLASSNAME::
506  _get_tensor3( size_t index_ ) const
507  {
508  const tensor3<I1, I2, I3, T>* tens3 = new tensor3<I1, I2, I3, T>( (void *)(_array + I1 * I2 * I3 * index_));
509  return *tens3;
510  }
511 
512 
513 
514  VMML_TEMPLATE_STRING
515  inline void
516  VMML_TEMPLATE_CLASSNAME::
517  get_tensor3( const size_t i4_, tensor3_t& t3_data_ ) const
518  {
519 #ifdef VMMLIB_SAFE_ACCESSORS
520  if ( i4_ >= I4 )
521  VMMLIB_ERROR( "get_tensor3() - index out of bounds.", VMMLIB_HERE );
522 #endif
523 
524  t3_data_ = _get_tensor3( i4_ );
525  }
526 
527 
528 
529  VMML_TEMPLATE_STRING
530  inline typename VMML_TEMPLATE_CLASSNAME::tensor3_t&
531  VMML_TEMPLATE_CLASSNAME::
532  get_tensor3( size_t i4_ )
533  {
534 #ifdef VMMLIB_SAFE_ACCESSORS
535  if ( i4_ >= I4 )
536  VMMLIB_ERROR( "get_tensor3() - index out of bounds.", VMMLIB_HERE );
537 #endif
538  return _get_tensor3( i4_ );
539  }
540 
541 
542  VMML_TEMPLATE_STRING
543  inline const typename VMML_TEMPLATE_CLASSNAME::tensor3_t&
544  VMML_TEMPLATE_CLASSNAME::
545  get_tensor3( size_t i4_ ) const
546  {
547 #ifdef VMMLIB_SAFE_ACCESSORS
548  if ( i4_ >= I4 )
549  VMMLIB_ERROR( "get_tensor3() - index out of bounds.", VMMLIB_HERE );
550 #endif
551  return _get_tensor3( i4_ );
552  }
553 
554 
555 
556  VMML_TEMPLATE_STRING
557  inline void
558  VMML_TEMPLATE_CLASSNAME::
559  set_tensor3( size_t i4_, const tensor3_t& t3_data_ )
560  {
561 #ifdef VMMLIB_SAFE_ACCESSORS
562  if ( i4_ >= I4 )
563  VMMLIB_ERROR( "set_tensor3() - index out of bounds.", VMMLIB_HERE );
564 #endif
565  memcpy(_array + I1*I2*I3*i4_, t3_data_.get_array_ptr(), I1*I2*I3*sizeof(T));
566  }
567 
568 
569 
570  //fill
571  VMML_TEMPLATE_STRING
572  void VMML_TEMPLATE_CLASSNAME::fill( T fillValue )
573  {
574  for( size_t i4 = 0; i4 < I4; ++i4 )
575  {
576  for( size_t i3 = 0; i3 < I3; ++i3 )
577  {
578  for( size_t i1 = 0; i1 < I1; ++i1 )
579  {
580  for( size_t i2 = 0; i2 < I2; ++i2 )
581  {
582  at(i1, i2, i3, i4) = fillValue;
583  }
584  }
585  }
586  }
587  //FIXME:
588  //_get_tensor3( i4 ).fill( fillValue );
589 
590  }
591 
592 
593  VMML_TEMPLATE_STRING
594  inline T&
595  VMML_TEMPLATE_CLASSNAME::at( size_t i1, size_t i2, size_t i3, size_t i4 )
596  {
597 #ifdef VMMLIB_SAFE_ACCESSORS
598  if ( i1 >= I1 || i2 >= I2 || i3 >= I3 || i4 >= I4 )
599  VMMLIB_ERROR( "at( i1, i2, i3, i4 ) - index out of bounds", VMMLIB_HERE );
600 #endif
601  return _array[ i4 * T3_SIZE + i3 * MATRIX_SIZE + i2 * ROWS + i1 ];
602 
603  }
604 
605 
606 
607  VMML_TEMPLATE_STRING
608  const inline T&
609  VMML_TEMPLATE_CLASSNAME::at( size_t i1, size_t i2, size_t i3, size_t i4 ) const
610  {
611 #ifdef VMMLIB_SAFE_ACCESSORS
612  if ( i1 >= I1 || i2 >= I2 || i3 >= I3 || i4 >= I4 )
613  VMMLIB_ERROR( "at( i1, i2, i3, i4 ) - index out of bounds", VMMLIB_HERE );
614 #endif
615  return _array[ i4 * T3_SIZE + i3 * MATRIX_SIZE + i2 * ROWS + i1 ];
616  }
617 
618 
619 
620  VMML_TEMPLATE_STRING
621  inline T&
622  VMML_TEMPLATE_CLASSNAME::operator()( size_t i1, size_t i2, size_t i3, size_t i4 )
623  {
624  return at( i1, i2, i3, i4 );
625  }
626 
627 
628 
629  VMML_TEMPLATE_STRING
630  const inline T&
631  VMML_TEMPLATE_CLASSNAME::operator()( size_t i1, size_t i2, size_t i3, size_t i4 ) const
632  {
633  return at( i1, i2, i3, i4 );
634  }
635 
636  /*VMML_TEMPLATE_STRING
637  void
638  VMML_TEMPLATE_CLASSNAME::fill( T fill_value )
639  {
640  for(size_t index = 0; index < SIZE; ++index)
641  {
642  _array[ index ] = fill_value;
643  }
644  }*/
645 
646  VMML_TEMPLATE_STRING
647  void
648  VMML_TEMPLATE_CLASSNAME::operator=( T fill_value )
649  {
650  fill( fill_value );
651  }
652 
653  VMML_TEMPLATE_STRING
654  void
655  VMML_TEMPLATE_CLASSNAME::zero()
656  {
657  fill( static_cast< T >( 0.0 ) );
658  }
659 
660  VMML_TEMPLATE_STRING
661  void
662  VMML_TEMPLATE_CLASSNAME::fill_increasing_values()
663  {
664  double fillValue = 0.0f;
665  for (size_t i4 = 0; i4 < I4; ++i4) {
666  for (size_t i3 = 0; i3 < I3; ++i3) {
667  for (size_t i1 = 0; i1 < I1; ++i1) {
668  for (size_t i2 = 0; i2 < I2; ++i2) {
669  at(i1, i2, i3, i4) = static_cast<T> (fillValue);
670  fillValue++;
671  }
672  }
673  }
674  }
675  }
676 
677  VMML_TEMPLATE_STRING
678  void
679  VMML_TEMPLATE_CLASSNAME::fill_random(int seed)
680  {
681  if ( seed >= 0 )
682  srand( seed );
683 
684  double fillValue = 0.0f;
685  for( size_t index = 0; index < SIZE; ++index )
686  {
687  fillValue = rand();
688  fillValue /= RAND_MAX;
689  fillValue *= (std::numeric_limits< T >::max)();
690  _array[ index ] = static_cast< T >( fillValue ) ;
691  }
692  }
693 
694  VMML_TEMPLATE_STRING
695  void
696  VMML_TEMPLATE_CLASSNAME::fill_random_signed(int seed)
697  {
698  if ( seed >= 0 )
699  srand( seed );
700 
701  double fillValue = 0.0f;
702  for( size_t index = 0; index < SIZE; ++index )
703  {
704 
705  fillValue = rand();
706  fillValue /= RAND_MAX;
707  fillValue *= (std::numeric_limits< T >::max)();
708  T fillValue2 = static_cast< T >(fillValue) % (std::numeric_limits< T >::max)();
709  fillValue2 -= (std::numeric_limits< T >::max)()/2;
710  _array[ index ] = fillValue2 ;
711  // test if ever > max/2 or < -max/2
712 
713  }
714  }
715 
716  VMML_TEMPLATE_STRING
717  const VMML_TEMPLATE_CLASSNAME&
718  VMML_TEMPLATE_CLASSNAME::operator=( const VMML_TEMPLATE_CLASSNAME& source_ )
719  {
720  if(this != &source_) // avoid self assignment
721  {
722  memcpy( _array, source_._array, I1 * I2 * I3 * I4 * sizeof( T ) );
723  }
724  return *this;
725  }
726 
727  VMML_TEMPLATE_STRING
728  template< typename TT >
729  void
730  VMML_TEMPLATE_CLASSNAME::cast_from( const tensor4< I1, I2, I3, I4, TT >& other )
731  {
732 #pragma omp parallel for
733  for(long tensor_idx = 0; tensor_idx < (long)I4; ++tensor_idx)
734  {
735 #pragma omp parallel for
736  for( long slice_idx = 0; slice_idx < (long)I3; ++slice_idx )
737  {
738 #pragma omp parallel for
739  for( long row_index = 0; row_index < (long)I1; ++row_index )
740  {
741 #pragma omp parallel for
742  for( long col_index = 0; col_index < (long)I2; ++col_index )
743  {
744  at( row_index, col_index, slice_idx, tensor_idx ) = static_cast< T >(other.at( row_index, col_index, slice_idx, tensor_idx ));
745  }
746  }
747  }
748  }
749 
750  }
751 
752  VMML_TEMPLATE_STRING
753  template< size_t J1, size_t J2, size_t J3, size_t J4, typename TT >
754  void
755  VMML_TEMPLATE_CLASSNAME::cast_from( const tensor4< J1, J2, J3, J4, TT >& other )
756  {
757  this->zero();
758  int j4 = I4 < J4? I4 : J4;
759  int j3 = I3 < J3? I3 : J3;
760  int j2 = I2 < J2? I2 : J2;
761  int j1 = I1 < J1? I1 : J1;
762 #pragma omp parallel for
763  for(long tensor_idx = 0; tensor_idx < (long)j4; ++tensor_idx)
764  {
765 #pragma omp parallel for
766  for( long slice_idx = 0; slice_idx < (long)j3; ++slice_idx )
767  {
768 #pragma omp parallel for
769  for( long row_index = 0; row_index < (long)j1; ++row_index )
770  {
771 #pragma omp parallel for
772  for( long col_index = 0; col_index < (long)j2; ++col_index )
773  {
774  at( row_index, col_index, slice_idx, tensor_idx ) = static_cast< T >(other.at( row_index, col_index, slice_idx, tensor_idx ));
775  }
776  }
777  }
778  }
779 
780  }
781 
782  VMML_TEMPLATE_STRING
783  template< typename TT >
784  void
785  VMML_TEMPLATE_CLASSNAME::float_t_to_uint_t( const tensor4< I1, I2, I3, I4, TT >& other )
786  {
787  if( sizeof(T) == 1 || sizeof(T) == 2) {
788  const TT* otherdata = other.get_array_ptr();
789  for( size_t index = 0;index < SIZE; ++index )
790  {
791  _array[index] = T( (std::min)( (std::max)(int(0), int( otherdata[index] + 0.5)), int((std::numeric_limits< T >::max)()) ));
792  }
793  }
794  else {
795  this->cast_from( other );
796  return;
797  }
798  }
799 
800  VMML_TEMPLATE_STRING
801  inline VMML_TEMPLATE_CLASSNAME
802  VMML_TEMPLATE_CLASSNAME::operator+( T scalar ) const
803  {
804  vmml::tensor4<I1, I2, I3, I4, T> result(*this);
805  result += scalar;
806  return result;
807 
808  }
809 
810  VMML_TEMPLATE_STRING
811  inline VMML_TEMPLATE_CLASSNAME
812  VMML_TEMPLATE_CLASSNAME::operator-( T scalar ) const
813  {
814  vmml::tensor4<I1, I2, I3, I4, T> result(*this);
815  result -= scalar;
816  return result;
817 
818  }
819 
820  VMML_TEMPLATE_STRING
821  void
822  VMML_TEMPLATE_CLASSNAME::operator+=( T scalar )
823  {
824  for(size_t index = 0; index < SIZE; ++index)
825  {
826  _array[index] += scalar;
827  }
828 
829  }
830 
831  VMML_TEMPLATE_STRING
832  void
833  VMML_TEMPLATE_CLASSNAME::operator-=( T scalar )
834  {
835  for(size_t index = 0; index < SIZE; ++index)
836  {
837  _array[index] -= scalar;
838  }
839 
840  }
841 
842  VMML_TEMPLATE_STRING
843  inline VMML_TEMPLATE_CLASSNAME
844  VMML_TEMPLATE_CLASSNAME::operator+( const tensor4& other ) const
845  {
846  vmml::tensor4<I1, I2, I3, I4, T> result(*this);
847  result += other;
848  return result;
849  }
850 
851 
852  VMML_TEMPLATE_STRING
853  inline VMML_TEMPLATE_CLASSNAME
854  VMML_TEMPLATE_CLASSNAME::operator-( const tensor4& other ) const
855  {
856  vmml::tensor4<I1, I2, I3, I4, T> result(*this);
857  result -= other;
858  return result;
859  }
860 
861  VMML_TEMPLATE_STRING
862  void
863  VMML_TEMPLATE_CLASSNAME::operator+=( const tensor4& other )
864  {
865  const T* dataptr = other.get_array_ptr();
866  for(size_t index = 0; index < SIZE; ++index)
867  {
868  _array[index] += dataptr[index];
869  }
870  }
871 
872  VMML_TEMPLATE_STRING
873  void
874  VMML_TEMPLATE_CLASSNAME::operator-=( const tensor4& other )
875  {
876  const T* dataptr = other.get_array_ptr();
877  for(size_t index = 0; index < SIZE; ++index)
878  {
879  _array[index] -= dataptr[index];
880  }
881  }
882 
883  VMML_TEMPLATE_STRING
884  template< size_t J1, size_t J2, size_t J3, size_t J4 >
885  typename enable_if< J1 < I1 && J2 < I2 && J3 < I3 && J4 < I4 >::type*
886  VMML_TEMPLATE_CLASSNAME::operator+=( const tensor4< J1, J2, J3, J4, T>& other )
887  {
888  for(size_t i4 = 0; i4 < J4; ++i4)
889  {
890  for( size_t i3 = 0; i3 < J3; ++i3 )
891  {
892  for( size_t i1 = 0; i1 < J1; ++i1 )
893  {
894  for( size_t i2 = 0; i2 < J2; ++i2 )
895  {
896  at( i1, i2, i3, i4 ) += other.at(i1, i2, i3, i4);
897  }
898  }
899  }
900  }
901  return 0;
902  }
903 
904 
905 
906  //
907  // tensor4-scalar operations / scaling
908  //
909  VMML_TEMPLATE_STRING
910  VMML_TEMPLATE_CLASSNAME
911  VMML_TEMPLATE_CLASSNAME::operator*( T scalar )
912  {
913  vmml::tensor4<I1, I2, I3, I4, T> result(*this);
914  result *= scalar;
915  return result;
916  }
917 
918  VMML_TEMPLATE_STRING
919  void
920  VMML_TEMPLATE_CLASSNAME::operator*=( T scalar )
921  {
922  for(size_t index = 0; index < SIZE; ++index)
923  {
924  _array[index] *= scalar;
925  }
926  }
927 
928  VMML_TEMPLATE_STRING
929  VMML_TEMPLATE_CLASSNAME
930  VMML_TEMPLATE_CLASSNAME::operator/( T scalar )
931  {
932  vmml::tensor4<I1, I2, I3, I4, T> result(*this);
933  result /= scalar;
934  return result;
935  }
936 
937  VMML_TEMPLATE_STRING
938  void
939  VMML_TEMPLATE_CLASSNAME::operator/=( T scalar )
940  {
941  for(size_t index = 0; index < SIZE; ++index)
942  {
943  _array[index] /= scalar;
944  }
945  }
946 
947 
948  VMML_TEMPLATE_STRING
949  inline VMML_TEMPLATE_CLASSNAME
950  VMML_TEMPLATE_CLASSNAME::operator-() const
951  {
952  return negate();
953  }
954 
955  VMML_TEMPLATE_STRING
956  VMML_TEMPLATE_CLASSNAME
957  VMML_TEMPLATE_CLASSNAME::negate() const
958  {
959  vmml::tensor4<I1, I2, I3, I4, T> result(*this);
960  T* dataptr = result.get_array_ptr();
961 
962  for(size_t index = 0; index < SIZE; ++index)
963  {
964  dataptr[index] = _array[index] * -1;
965  }
966  return result;
967 
968  }
969 
970  VMML_TEMPLATE_STRING
971  typename VMML_TEMPLATE_CLASSNAME::tensor3_t&
972  VMML_TEMPLATE_CLASSNAME::average_I4(tensor3_t& t3) const
973  {
974 
975  for(size_t i3 = 0; i3 < I3; ++i3)
976  {
977  for( size_t i1 = 0; i1 < I1; ++i1 )
978  {
979  for( size_t i2 = 0; i2 < I2; ++i2 )
980  {
981  T sum = 0;
982  for( size_t i4 = 0; i4 < I4; ++i4 )
983  {
984  sum += at(i1,i2,i3,i4);
985  }
986  t3.at(i1,i2,i3) = sum/I4;
987  }
988  }
989  }
990  return t3;
991 
992  }
993 
994  VMML_TEMPLATE_STRING
995  size_t
996  VMML_TEMPLATE_CLASSNAME::get_array_size_in_bytes()
997  {
998  return (sizeof(T) * SIZE);
999  }
1000 
1001 
1002 
1003  // WARNING: dangerous. Use before destruction if you want to prevent
1004  // a delete call for the assigned T* _array in the destructor.
1005  VMML_TEMPLATE_STRING
1006  void
1007  VMML_TEMPLATE_CLASSNAME::clear_array_pointer()
1008  {
1009  _array = 0;
1010  }
1011 
1012 // VMML_TEMPLATE_STRING
1013 // void
1014 // VMML_TEMPLATE_CLASSNAME::mode1_unfolding_fwd(mode1_unfolding_type& unfolding) const {
1015 // typedef matrix< I1, I2, T > slice_type;
1016 // slice_type* slice = new slice_type();
1017 // for (size_t l = 0; l < I4; ++l) {
1018 // tensor3_t t3;
1019 // get_tensor3(l,t3);
1020 // for (size_t k = 0; k < I3; ++k) {
1021 // t3.get_frontal_slice_fwd(k, *slice);
1022 // for (size_t j = 0; j < I2; ++j) {
1023 // unfolding.set_column(l*I2*I3 + k*I2 + j, slice->get_column(j));
1024 // }
1025 // }
1026 // }
1027 // delete slice;
1028 // }
1029 //
1030 // VMML_TEMPLATE_STRING
1031 // void
1032 // VMML_TEMPLATE_CLASSNAME::mode2_unfolding_fwd(mode2_unfolding_type& unfolding) const {
1033 // typedef matrix< I2, I3, T > slice_type;
1034 // slice_type* slice = new slice_type();
1035 // for (size_t i = 0; i < I1; ++i) {
1036 // for (size_t l = 0; l < I4; ++l) {
1037 // tensor3_t t3;
1038 // get_tensor3(l,t3);
1039 // t3.get_horizontal_slice_fwd(i, *slice);
1040 // for (size_t k = 0; k < I3; ++k) {
1041 // unfolding.set_column(i*I3*I4 + l*I3 + k, slice->get_column(k));
1042 // }
1043 // }
1044 // }
1045 // delete slice;
1046 // }
1047 //
1048 // VMML_TEMPLATE_STRING
1049 // void
1050 // VMML_TEMPLATE_CLASSNAME::mode3_unfolding_fwd(mode3_unfolding_type& unfolding) const {
1051 // typedef matrix< I2, I3, T > slice_type;
1052 // slice_type* slice = new slice_type();
1053 // for (size_t j = 1; j < I2; ++j) {
1054 // for (size_t i = 0; i < I1; ++i) {
1055 // for (size_t l = 0; l < I4; ++l) {
1056 // tensor3_t t3;
1057 // get_tensor3(l,t3);
1058 // t3.get_horizontal_slice_fwd(i, *slice);
1059 // unfolding.set_column(j*I4*I1 + i*I4 + l, slice->get_column(j));
1060 // }
1061 // }
1062 // }
1063 // delete slice;
1064 // }
1065 //
1066 // VMML_TEMPLATE_STRING
1067 // void
1068 // VMML_TEMPLATE_CLASSNAME::mode4_unfolding_fwd(mode4_unfolding_type& unfolding) const {
1069 // for (size_t k = 1; k < I3; ++k) {
1070 // for (size_t j = 0; j < I2; ++j) {
1071 // for (size_t i = 0; i < I1; ++i) {
1072 // for (size_t l = 0; l < I4; ++l) {
1073 // unfolding.at(l,k*I1*I2 + j*I1 + i) = at(i,j,k,l);
1074 // }
1075 // }
1076 // }
1077 // }
1078 // }
1079 
1080  VMML_TEMPLATE_STRING
1081  void
1082  VMML_TEMPLATE_CLASSNAME::mode1_unfolding_fwd(mode1_unfolding_type& unfolding) const {
1083  typedef matrix< I1, I2, T > slice_type;
1084  slice_type* slice = new slice_type();
1085  for (size_t l = 0; l < I4; ++l) {
1086  tensor3_t t3;
1087  get_tensor3(l,t3);
1088  for (size_t k = 0; k < I3; ++k) {
1089  t3.get_frontal_slice_fwd(k, *slice);
1090  for (size_t j = 0; j < I2; ++j) {
1091  unfolding.set_column(l*I2*I3 + k*I2 + j, slice->get_column(j));
1092  }
1093  }
1094  }
1095  delete slice;
1096  }
1097 
1098  VMML_TEMPLATE_STRING
1099  void
1100  VMML_TEMPLATE_CLASSNAME::mode2_unfolding_fwd(mode2_unfolding_type& unfolding) const {
1101  typedef matrix< I2, I3, T > slice_type;
1102  slice_type* slice = new slice_type();
1103  for (size_t i = 0; i < I1; ++i) {
1104  for (size_t l = 0; l < I4; ++l) {
1105  tensor3_t t3;
1106  get_tensor3(l,t3);
1107  t3.get_horizontal_slice_fwd(i, *slice);
1108  for (size_t k = 0; k < I3; ++k) {
1109  unfolding.set_column(i*I3*I4 + l*I3 + k, slice->get_column(k));
1110  }
1111  }
1112  }
1113  delete slice;
1114  }
1115 
1116  VMML_TEMPLATE_STRING
1117  void
1118  VMML_TEMPLATE_CLASSNAME::mode3_unfolding_fwd(mode3_unfolding_type& unfolding) const {
1119  typedef matrix< I2, I3, T > slice_type;
1120  slice_type* slice = new slice_type();
1121  for (size_t j = 0; j < I2; ++j) {
1122  for (size_t i = 0; i < I1; ++i) {
1123  for (size_t l = 0; l < I4; ++l) {
1124  tensor3_t t3;
1125  get_tensor3(l,t3);
1126  t3.get_horizontal_slice_fwd(i, *slice);
1127  unfolding.set_column(j*I4*I1 + i*I4 + l, slice->get_row(j));
1128  }
1129  }
1130  }
1131  delete slice;
1132  }
1133 
1134  VMML_TEMPLATE_STRING
1135  void
1136  VMML_TEMPLATE_CLASSNAME::mode4_unfolding_fwd(mode4_unfolding_type& unfolding) const {
1137  for (size_t k = 0; k < I3; ++k) {
1138  for (size_t j = 0; j < I2; ++j) {
1139  for (size_t i = 0; i < I1; ++i) {
1140  for (size_t l = 0; l < I4; ++l) {
1141  unfolding.at(l,k*I1*I2 + j*I1 + i) = at(i,j,k,l);
1142  }
1143  }
1144  }
1145  }
1146  }
1147 
1148  VMML_TEMPLATE_STRING
1149  double
1150  VMML_TEMPLATE_CLASSNAME::frobenius_norm() const {
1151  double f_norm = 0.0;
1152  for (long i3 = 0; i3 < long(I3); ++i3) {
1153  for (long i4 = 0; i4 < long(I4); ++i4) {
1154  for (long i1 = 0; i1 < long(I1); ++i1) {
1155  long i2 = 0;
1156  for (i2 = 0; i2 < long(I2); ++i2) {
1157  f_norm += at(i1, i2, i3, i4) * at(i1, i2, i3, i4);
1158  }
1159  }
1160  }
1161  }
1162  return sqrt(f_norm);
1163  }
1164 
1165  VMML_TEMPLATE_STRING
1166  double
1167  VMML_TEMPLATE_CLASSNAME::frobenius_norm(const tensor4< I1, I2, I3, I4, T>& other_) const {
1168  double f_norm = 0.0;
1169  T abs_diff = 0;
1170  T *it = _array, *it_end = _array + I1*I2*I3*I4;
1171  T *it_other = other_._array;
1172  for (; it != it_end; ++it, ++it_other) {
1173  abs_diff = fabs(*it - *it_other);
1174  f_norm += abs_diff * abs_diff;
1175  }
1176 
1177  return sqrt(f_norm);
1178  }
1179 
1180  VMML_TEMPLATE_STRING
1181  template< typename TT >
1182  void
1183  VMML_TEMPLATE_CLASSNAME::quantize_log(tensor4< I1, I2, I3, I4, TT >& quantized_, tensor4< I1, I2, I3, I4, char >& signs_, T& min_value_, T& max_value_, const TT& tt_range_) const {
1184  double max_tt_range = double(tt_range_);
1185  double min_tt_range = 0;
1186 
1187  min_value_ = 0;
1188  max_value_ = get_abs_max();
1189  double t_range = max_value_ - min_value_;
1190 
1191  T* it = _array;
1192  TT* it_quant = quantized_.get_array_ptr();
1193  char* it_sign = signs_.get_array_ptr();
1194  for (; it != _array+I1*I2*I3*I4; ++it, ++it_quant, ++it_sign) {
1195  T value = fabs(*it);
1196  *it_sign = ((*it) < 0.f) ? 0 : 1;
1197  T quant_value = 0;
1198  if (std::numeric_limits<TT>::is_signed) {
1199  quant_value = log2(1 + value) / log2(1 + t_range) * tt_range_;
1200  *it_quant = TT((std::min)((std::max)(min_tt_range, double(quant_value + 0.5)), max_tt_range));
1201  } else {
1202  quant_value = log2(1 + (value - min_value_)) / log2(1 + t_range) * tt_range_;
1203  *it_quant = TT((std::min)((std::max)(min_tt_range, double(quant_value + 0.5)), max_tt_range));
1204  }
1205  }
1206  }
1207 
1208  VMML_TEMPLATE_STRING
1209  double
1210  VMML_TEMPLATE_CLASSNAME::mean() const {
1211  double val = 0;
1212  T* it = _array, *it_end = _array + I1*I2*I3*I4;
1213  for (; it != it_end; ++it) {
1214  val += double(abs(*it));
1215  }
1216 
1217  return ( val / size());
1218  }
1219 
1220  VMML_TEMPLATE_STRING
1221  void
1222  VMML_TEMPLATE_CLASSNAME::mean(T& mean_) const {
1223  mean_ = static_cast<T> (mean());
1224  }
1225 
1226  VMML_TEMPLATE_STRING
1227  double
1228  VMML_TEMPLATE_CLASSNAME::variance() const {
1229  double val = 0.0;
1230  double sum_val = 0.0;
1231  double mean_val = mean();
1232  T* it = _array, *it_end = _array+I1*I2*I3*I4;
1233  for (; it != it_end; ++it) {
1234  val = double(*it) - mean_val;
1235  val *= val;
1236  sum_val += val;
1237  }
1238 
1239  return double(sum_val / (size() - 1));
1240  }
1241 
1242  VMML_TEMPLATE_STRING
1243  double
1244  VMML_TEMPLATE_CLASSNAME::stdev() const {
1245  return sqrt(variance());
1246  }
1247 
1248  VMML_TEMPLATE_STRING
1249  T
1250  VMML_TEMPLATE_CLASSNAME::get_min() const {
1251  T tensor4_min = static_cast<T> ((std::numeric_limits<T>::max)());
1252 
1253  T *it = _array, *it_end = _array + I1*I2*I3*I4;
1254  for (; it != it_end; ++it) {
1255  if (*it < tensor4_min) {
1256  tensor4_min = *it;
1257  }
1258  }
1259  return tensor4_min;
1260  }
1261 
1262  VMML_TEMPLATE_STRING
1263  T
1264  VMML_TEMPLATE_CLASSNAME::get_max() const {
1265  T tensor4_max = static_cast<T> (0);
1266 
1267  T *it = _array, *it_end = _array + I1*I2*I3*I4;
1268  for (; it != it_end; ++it) {
1269  if (*it > tensor4_max) {
1270  tensor4_max = *it;
1271  }
1272  }
1273  return tensor4_max;
1274  }
1275 
1276  VMML_TEMPLATE_STRING
1277  T
1278  VMML_TEMPLATE_CLASSNAME::get_abs_min() const {
1279  T tensor4_min = static_cast<T> ((std::numeric_limits<T>::max)());
1280 
1281  T *it = _array, *it_end = _array + I1*I2*I3*I4;
1282  for (; it != it_end; ++it) {
1283  if (fabs(*it) < fabs(tensor4_min)) {
1284  tensor4_min = fabs(*it);
1285  }
1286  }
1287  return tensor4_min;
1288  }
1289 
1290  VMML_TEMPLATE_STRING
1291  T
1292  VMML_TEMPLATE_CLASSNAME::get_abs_max() const {
1293  T tensor4_max = static_cast<T> (0);
1294 
1295  T *it = _array, *it_end = _array + I1*I2*I3*I4;
1296  for (; it != it_end; ++it) {
1297  if (fabs(*it) > fabs(tensor4_max)) {
1298  tensor4_max = fabs(*it);
1299  }
1300  }
1301  return tensor4_max;
1302  }
1303 
1304 #undef VMML_TEMPLATE_STRING
1305 #undef VMML_TEMPLATE_CLASSNAME
1306 
1307  } // namespace vmml
1308 
1309 #endif