Lib/C2/lm/vmatrix.h

   1 #pragma once
   2
   3 #include <vector>
   4 #include <exception>
   5
   6 #include <cassert>
   7
   8
   9
  10 namespace lm
  11 {
  12
  13
  14 class invalid_product_exception : std::exception
  15 {};
  16
  17
  18 //! 可変サイズ行列.
  19 /*
  20  size_x == 3 , size_y == 2 とすると
  21  [ (0,0) , (1,0) , (2,0) ]
  22  [ (0,1) , (1,1) , (2,1) ]
  23
  24  通しインデックスは
  25  [ 0 , 1 , 2 ]
  26  [ 3 , 4 , 5 ]
  27 */
  28 template<typename T>
  29 class matrix
  30 {
  31 public:
  32         matrix(void);
  33         matrix(size_t num_x, size_t num_y);
  34         matrix(size_t num_x, size_t num_y, const T& initial_value);
  35         matrix(const matrix<T>& m);
  36
  37         T* v(void);
  38         const T* v(void) const;
  39
  40         T& at(size_t idx_v);
  41         const T& at(size_t idx_v) const;
  42         T& operator()(size_t idx_x, size_t idx_y);
  43         const T& operator()(size_t idx_x, size_t idx_y) const;
  44         T& operator[](size_t idx_v);
  45         const T& operator[](size_t idx_v) const;
  46
  47         T& at(size_t idx_x, size_t idx_y);
  48         const T& at(size_t idx_x, size_t idx_y) const;
  49         T& operator()(size_t idx_v);
  50         const T& operator()(size_t idx_v) const;
  51
  52         void resize(size_t num_x, size_t num_y);
  53         void resize(size_t num_x, size_t num_y, const T& initial_value);
  54
  55         void clear(void);
  56
  57         size_t size(void) const;
  58         size_t size_x(void) const;
  59         size_t size_y(void) const;
  60
  61         bool empty(void) const;
  62
  63         void set(const matrix<T>& m);
  64         matrix<T>& operator=(const matrix<T>& m);
  65
  66         bool equals(const matrix<T>& m);
  67         bool operator==(const matrix<T>& m);
  68         bool operator!=(const matrix<T>& m);
  69
  70         typedef typename std::vector<T>::iterator               iterator;
  71         typedef typename std::vector<T>::const_iterator         const_iterator;
  72         typedef typename std::vector<T>::reverse_iterator       reverse_iterator;
  73         typedef typename std::vector<T>::const_reverse_iterator const_reverse_iterator;
  74         iterator               begin  (void);
  75         iterator               end    (void);
  76         reverse_iterator       rbegin (void);
  77         reverse_iterator       rend   (void);
  78         const_iterator         begin  (void) const;
  79         const_iterator         end    (void) const;
  80         const_reverse_iterator rbegin (void) const;
  81         const_reverse_iterator rend   (void) const;
  82
  83         void swap(matrix<T>& m);
  84
  85 protected:
  86         size_t m_size_x;
  87         size_t m_size_y;
  88         std::vector<T> buf;
  89 };
  90
  91
  92 // global method
  93
  94 template<typename T> inline
  95 std::vector<T> dot(const matrix<T>& _m, const std::vector<T>& _v);
  96 template<typename T> inline
  97 std::vector<T> dot(const std::vector<T>& _v, const matrix<T>& _m);
  98
  99 template<typename T> inline
 100 std::vector<T> operator*(const matrix<T>& _m, const std::vector<T>& _v);
 101 template<typename T> inline
 102 std::vector<T> operator*(const std::vector<T>& _v, const matrix<T>& _m);
 103
 104 template<typename T> inline
 105 matrix<T> dot(const matrix<T>& _m0, const matrix<T>& _m1);
 106 template<typename T> inline
 107 matrix<T> operator*(const matrix<T>& _m0, const matrix<T>& _m1);
 108
 109
 110 // implement
 111
 112 template<typename T> inline
 113 matrix<T>::matrix(void)
 114         : m_size_x( 0 )
 115         , m_size_y( 0 )
 116 {}
 117
 118 template<typename T> inline
 119 matrix<T>::matrix(size_t num_x, size_t num_y)
 120         : m_size_x( num_x )
 121         , m_size_y( num_y )
 122 {
 123         buf.resize( m_size_x * m_size_y );
 124 }
 125
 126 template<typename T> inline
 127 matrix<T>::matrix(size_t num_x, size_t num_y, const T& initial_value)
 128         : m_size_x( num_x )
 129         , m_size_y( num_y )
 130 {
 131         buf.resize( m_size_x * m_size_y );
 132         std::fill( buf.begin() , buf.end() , initial_value );
 133 }
 134
 135 template<typename T> inline
 136 matrix<T>::matrix(const matrix<T>& m)
 137 {
 138         buf = m.buf;
 139         m_size_x = m.size_x();
 140         m_size_y = m.size_y();
 141 }
 142
 143
 144 template<typename T> inline
 145 T* matrix<T>::v()
 146 {
 147         if( buf.empty() ) return NULL    ;
 148         else              return &buf[0] ;
 149 }
 150 template<typename T> inline
 151 const T* matrix<T>::v() const
 152 {
 153         if( buf.empty() ) return NULL    ;
 154         else              return &buf[0] ;
 155 }
 156
 157
 158 template<typename T> inline
 159 T& matrix<T>::at(size_t idx_v)
 160 {
 161         return buf[idx_v];
 162 }
 163 template<typename T> inline
 164 const T& matrix<T>::at(size_t idx_v) const
 165 {
 166         return buf[idx_v];
 167 }
 168 template<typename T> inline
 169 T& matrix<T>::operator()(size_t idx_v)
 170 {
 171         return at(idx_v);
 172 }
 173 template<typename T> inline
 174 const T& matrix<T>::operator()(size_t idx_v) const
 175 {
 176         return at(idx_v);
 177 }
 178 template<typename T> inline
 179 T& matrix<T>::operator[](size_t idx_v)
 180 {
 181         return at( idx_v );
 182 }
 183 template<typename T> inline
 184 const T& matrix<T>::operator[](size_t idx_v) const
 185 {
 186         return at( idx_v );
 187 }
 188
 189
 190 template<typename T> inline
 191 T& matrix<T>::at(size_t idx_x, size_t idx_y)
 192 {
 193         // Matrix subscript out of range
 194         assert( idx_x < m_size_x && idx_y < m_size_y );
 195
 196         return buf[ idx_x + m_size_x * idx_y ];
 197 }
 198 template<typename T> inline
 199 const T& matrix<T>::at(size_t idx_x, size_t idx_y) const
 200 {
 201         // Matrix subscript out of range
 202         assert( idx_x < m_size_x && idx_y < m_size_y );
 203
 204         return buf[ idx_x + m_size_x * idx_y ];
 205 }
 206 template<typename T> inline
 207 T& matrix<T>::operator()(size_t idx_x,size_t idx_y)
 208 {
 209         return at( idx_x , idx_y );
 210 }
 211 template<typename T> inline
 212 const T& matrix<T>::operator()(size_t idx_x, size_t idx_y) const
 213 {
 214         return at( idx_x , idx_y );
 215 }
 216
 217
 218 template<typename T> inline
 219 void matrix<T>::resize(size_t num_x, size_t num_y)
 220 {
 221         m_size_x = num_x ;
 222         m_size_y = num_y ;
 223         buf.resize( m_size_x * m_size_y );
 224 }
 225 template<typename T> inline
 226 void matrix<T>::resize(size_t num_x, size_t num_y, const T& initial_value)
 227 {
 228         m_size_x = num_x ;
 229         m_size_y = num_y ;
 230         buf.resize( m_size_x * m_size_y );
 231         std::fill( buf.begin() , buf.end() , initial_value );
 232 }
 233
 234
 235 template<typename T> inline
 236 void matrix<T>::clear(void)
 237 {
 238         m_size_x = 0;
 239         m_size_y = 0;
 240         buf.clear();
 241 }
 242
 243
 244 template<typename T> inline
 245 size_t matrix<T>::size(void) const
 246 {
 247         return buf.size();
 248 }
 249 template<typename T> inline
 250 size_t matrix<T>::size_x(void) const
 251 {
 252         return m_size_x;
 253 }
 254 template<typename T> inline
 255 size_t matrix<T>::size_y(void) const
 256 {
 257         return m_size_y;
 258 }
 259
 260 template<typename T> inline
 261 bool matrix<T>::empty(void) const
 262 {
 263         return buf.empty();
 264 }
 265
 266
 267 template<typename T> inline
 268 void matrix<T>::set(const matrix<T>& m)
 269 {
 270         buf = m.buf;
 271         m_size_x = m.size_x() ;   m_size_y = m.size_y() ;
 272 }
 273 template<typename T> inline
 274 matrix<T>& matrix<T>::operator=(const matrix<T>& m)
 275 {
 276         set( m ) ;
 277         return (*this);
 278 }
 279
 280
 281 template<typename T> inline
 282 bool matrix<T>::equals(const matrix<T>& m)
 283 {
 284         if( m_size_x != m.m_size_x || m_size_y != m.m_size_y ) return false;
 285         for(size_t i=0;i<m.size();i++) if( (*this)(i) != m(i) ) return false;
 286         return true;
 287 }
 288 template<typename T> inline
 289 bool matrix<T>::operator==(const matrix<T>& m)
 290 {
 291         return equals( m );
 292 }
 293 template<typename T> inline
 294 bool matrix<T>::operator!=(const matrix<T>& m)
 295 {
 296         return !equals( m );
 297 }
 298
 299
 300 template<typename T> inline
 301 typename matrix<T>::iterator matrix<T>::begin(void)
 302 {
 303         return buf.begin();
 304 }
 305 template<typename T> inline
 306 typename matrix<T>::iterator matrix<T>::end(void)
 307 {
 308         return buf.end();
 309 }
 310 template<typename T> inline
 311 typename matrix<T>::reverse_iterator matrix<T>::rbegin(void)
 312 {
 313         return buf.rbegin();
 314 }
 315 template<typename T> inline
 316 typename matrix<T>::reverse_iterator matrix<T>::rend(void)
 317 {
 318         return buf.rend();
 319 }
 320 template<typename T> inline
 321 typename matrix<T>::const_iterator matrix<T>::begin(void) const
 322 {
 323         return buf.begin();
 324 }
 325 template<typename T> inline
 326 typename matrix<T>::const_iterator matrix<T>::end(void) const
 327 {
 328         return buf.end();
 329 }
 330 template<typename T> inline
 331 typename matrix<T>::const_reverse_iterator matrix<T>::rbegin(void) const
 332 {
 333         return buf.rbegin();
 334 }
 335 template<typename T> inline
 336 typename matrix<T>::const_reverse_iterator matrix<T>::rend(void) const
 337 {
 338         return buf.rend();
 339 }
 340
 341 template<typename T> inline
 342 void matrix<T>::swap(matrix<T>& m)
 343 {
 344         std::swap(m_size_x, m.m_size_x);
 345         std::swap(m_size_y, m.m_size_y);
 346         buf.swap(m.buf);
 347 }
 348
 349
 350 template<typename T> inline
 351 std::vector<T> dot(const matrix<T>& _m, const std::vector<T>& _v)
 352 {
 353         return _m * _v;
 354 }
 355
 356 template<typename T> inline
 357 std::vector<T> dot(const std::vector<T>& _v, const matrix<T>& _m)
 358 {
 359         return _v * _m;
 360 }
 361
 362 //! 行列とベクトルの積
 363 //! m.size_x() == 3 , m.size_y() == 5 の場合
 364 //! [r0]   [m00 m10 m20]
 365 //! [r1]   [m01 m11 m21]   [v0]
 366 //! [r2] = [m02 m12 m22] * [v1]
 367 //! [r3]   [m03 m13 m23]   [v2]
 368 //! [r4]   [m04 m14 m24]
 369 //!
 370 //! ※ mij = m( i , j ) = m.at( i , j )
 371 template<typename T> inline
 372 std::vector<T> operator*(const matrix<T>& _m, const std::vector<T>& _v)
 373 {
 374         assert( _m.size_x() == _v.size() );
 375         if( _m.size_x() != _v.size() )
 376                 throw invalid_product_exception();
 377
 378         std::vector<T> ret( _m.size_y() );
 379         for( size_t i = 0 ; i < _m.size_y() ; ++i )
 380         {
 381                 T& val = ret[i];
 382
 383                 for( size_t j = 0 ; j < _m.size_x() ; ++j )
 384                 {
 385                         val += _m( j , i ) * _v[j];
 386                 }
 387         }
 388
 389         return ret;
 390 }
 391
 392 //! 行列とベクトルの積
 393 //! m.size_x() == 5 , m.size_y() == 3 の場合
 394 //!                                 [m00 m10 m20 m30 m40]
 395 //! [r0 r1 r2 r3 r4] = [v0 v1 v2] * [m01 m11 m21 m31 m41]
 396 //!                                 [m02 m12 m22 m32 m42]
 397 //!
 398 //! ※ mij = m( i , j ) = m.at( i , j )
 399 template<typename T> inline
 400 std::vector<T> operator*(const std::vector<T>& _v, const matrix<T>& _m)
 401 {
 402         assert( _m.size_y() == _v.size() );
 403         if( _m.size_y() != _v.size() )
 404                 throw invalid_product_exception();
 405
 406         std::vector<T> ret( _m.size_x() );
 407         for( size_t i = 0 ; i < _m.size_x() ; ++i )
 408         {
 409                 T& val = ret[i];
 410
 411                 for( size_t j = 0 ; j < _m.size_y() ; ++j )
 412                 {
 413                         val += _v[j] * _m( i , j );
 414                 }
 415         }
 416
 417         return ret;
 418 }
 419
 420 template<typename T> inline
 421 matrix<T> dot(const matrix<T>& _m0, const matrix<T>& _m1)
 422 {
 423         return _m0 * _m1;
 424 }
 425
 426 //! 行列同士の積
 427 //! A = M * N =>
 428 //! [a00 a10]   [m00 m10 m 20]
 429 //! [a01 a11]   [m01 m11 m 21]   [n00 n10]
 430 //! [a02 a12] = [m02 m12 m 22] * [n01 n11]
 431 //! [a03 a13]   [m03 m13 m 23]   [n02 n12]
 432 //! [a04 a14]   [m04 m14 m 24]
 433 template<typename T> inline
 434 matrix<T> operator*(const matrix<T>& _m0, const matrix<T>& _m1)
 435 {
 436         assert( _m0.size_x() == _m1.size_y() );
 437         if( _m0.size_x() != _m1.size_y() )
 438                 throw invalid_product_exception();
 439
 440         matrix<T> m( _m1.size_x() , _m0.size_y() );
 441         for( size_t i = 0 ; i < _m1.size_x() ; ++i )
 442         {
 443                 for( size_t j = 0 ; j < _m0.size_y() ; ++j )
 444                 {
 445                         T& val = m( i , j );
 446
 447                         for( size_t k = 1 ; k < _m0.size_x() ; ++k )
 448                         {
 449                                 val += _m0( k , j ) * _m1( i , k );
 450                         }
 451                 }
 452         }
 453
 454         return m;
 455 }
 456
 457
 458 }