src/2geom/numeric/matrix.h

   1 /*\r
   2  * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;\r
   3  * "views" mimic the semantic of C++ references: any operation performed\r
   4  * on a "view" is actually performed on the "viewed object"\r
   5  *\r
   6  * Authors:\r
   7  *              Marco Cecchetti <mrcekets at gmail.com>\r
   8  *\r
   9  * Copyright 2008  authors\r
  10  *\r
  11  * This library is free software; you can redistribute it and/or\r
  12  * modify it either under the terms of the GNU Lesser General Public\r
  13  * License version 2.1 as published by the Free Software Foundation\r
  14  * (the "LGPL") or, at your option, under the terms of the Mozilla\r
  15  * Public License Version 1.1 (the "MPL"). If you do not alter this\r
  16  * notice, a recipient may use your version of this file under either\r
  17  * the MPL or the LGPL.\r
  18  *\r
  19  * You should have received a copy of the LGPL along with this library\r
  20  * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
  21  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
  22  * You should have received a copy of the MPL along with this library\r
  23  * in the file COPYING-MPL-1.1\r
  24  *\r
  25  * The contents of this file are subject to the Mozilla Public License\r
  26  * Version 1.1 (the "License"); you may not use this file except in\r
  27  * compliance with the License. You may obtain a copy of the License at\r
  28  * http://www.mozilla.org/MPL/\r
  29  *\r
  30  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
  31  * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
  32  * the specific language governing rights and limitations.\r
  33  */\r
  34 \r
  35 \r
  36 \r
  37 \r
  38 #ifndef _NL_MATRIX_H_\r
  39 #define _NL_MATRIX_H_\r
  40 \r
  41 #include <2geom/numeric/vector.h>\r
  42 \r
  43 #include <cassert>\r
  44 #include <utility>    // for std::pair\r
  45 #include <algorithm>  // for std::swap\r
  46 #include <sstream>\r
  47 #include <string>\r
  48 \r
  49 #include <gsl/gsl_matrix.h>\r
  50 #include <gsl/gsl_linalg.h>\r
  51 \r
  52 \r
  53 namespace Geom { namespace NL {\r
  54 \r
  55 namespace detail\r
  56 {\r
  57 \r
  58 class BaseMatrixImpl\r
  59 {\r
  60   public:\r
  61         virtual ~BaseMatrixImpl()\r
  62         {\r
  63         }\r
  64 \r
  65         ConstVectorView row_const_view(size_t i) const\r
  66         {\r
  67                 return ConstVectorView(gsl_matrix_const_row(m_matrix, i));\r
  68         }\r
  69 \r
  70         ConstVectorView column_const_view(size_t i) const\r
  71         {\r
  72                 return ConstVectorView(gsl_matrix_const_column(m_matrix, i));\r
  73         }\r
  74 \r
  75         const double & operator() (size_t i, size_t j) const\r
  76         {\r
  77                 return *gsl_matrix_const_ptr(m_matrix, i, j);\r
  78         }\r
  79 \r
  80         const gsl_matrix* get_gsl_matrix() const\r
  81         {\r
  82                 return m_matrix;\r
  83         }\r
  84 \r
  85         bool is_zero() const\r
  86         {\r
  87                 return gsl_matrix_isnull(m_matrix);\r
  88         }\r
  89 \r
  90         bool is_positive() const\r
  91         {\r
  92                 return gsl_matrix_ispos(m_matrix);\r
  93         }\r
  94 \r
  95         bool is_negative() const\r
  96         {\r
  97                 return gsl_matrix_isneg(m_matrix);\r
  98         }\r
  99 \r
 100         bool is_non_negative() const\r
 101         {\r
 102                 for ( unsigned int i = 0; i < rows(); ++i )\r
 103                 {\r
 104                         for ( unsigned int j = 0; j < columns(); ++j )\r
 105                         {\r
 106                                 if ( (*this)(i,j) < 0 ) return false;\r
 107                         }\r
 108                 }\r
 109                 return true;\r
 110         }\r
 111 \r
 112         double max() const\r
 113         {\r
 114                 return gsl_matrix_max(m_matrix);\r
 115         }\r
 116 \r
 117         double min() const\r
 118         {\r
 119                 return gsl_matrix_min(m_matrix);\r
 120         }\r
 121 \r
 122         std::pair<size_t, size_t>\r
 123         max_index() const\r
 124         {\r
 125                 std::pair<size_t, size_t> indices;\r
 126                 gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));\r
 127                 return indices;\r
 128         }\r
 129 \r
 130         std::pair<size_t, size_t>\r
 131         min_index() const\r
 132         {\r
 133                 std::pair<size_t, size_t> indices;\r
 134                 gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));\r
 135                 return indices;\r
 136         }\r
 137 \r
 138         size_t rows() const\r
 139         {\r
 140                 return m_rows;\r
 141         }\r
 142 \r
 143         size_t columns() const\r
 144         {\r
 145                 return m_columns;\r
 146         }\r
 147 \r
 148         std::string str() const;\r
 149 \r
 150   protected:\r
 151         size_t m_rows, m_columns;\r
 152         gsl_matrix* m_matrix;\r
 153 \r
 154 };  // end class BaseMatrixImpl\r
 155 \r
 156 \r
 157 inline\r
 158 bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)\r
 159 {\r
 160         if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;\r
 161 \r
 162         for (size_t i = 0; i < m1.rows(); ++i)\r
 163                 for (size_t j = 0; j < m1.columns(); ++j)\r
 164                         if (m1(i,j) != m2(i,j)) return false;\r
 165 \r
 166         return true;\r
 167 }\r
 168 \r
 169 template< class charT >\r
 170 inline\r
 171 std::basic_ostream<charT> &\r
 172 operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix)\r
 173 {\r
 174         if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;\r
 175 \r
 176         os << "[[" << _matrix(0,0);\r
 177         for (size_t j = 1; j < _matrix.columns(); ++j)\r
 178         {\r
 179                 os << ", " << _matrix(0,j);\r
 180         }\r
 181         os << "]";\r
 182 \r
 183         for (size_t i = 1; i < _matrix.rows(); ++i)\r
 184         {\r
 185                 os << ", [" << _matrix(i,0);\r
 186                 for (size_t j = 1; j < _matrix.columns(); ++j)\r
 187                 {\r
 188                         os << ", " << _matrix(i,j);\r
 189                 }\r
 190                 os << "]";\r
 191         }\r
 192         os << "]";\r
 193         return os;\r
 194 }\r
 195 \r
 196 inline\r
 197 std::string BaseMatrixImpl::str() const\r
 198 {\r
 199         std::ostringstream oss;\r
 200         oss << (*this);\r
 201         return oss.str();\r
 202 }\r
 203 \r
 204 \r
 205 class MatrixImpl : public BaseMatrixImpl\r
 206 {\r
 207   public:\r
 208 \r
 209         typedef BaseMatrixImpl base_type;\r
 210 \r
 211         void set_all( double x )\r
 212         {\r
 213                 gsl_matrix_set_all(m_matrix, x);\r
 214         }\r
 215 \r
 216         void set_identity()\r
 217         {\r
 218                 gsl_matrix_set_identity(m_matrix);\r
 219         }\r
 220 \r
 221         using base_type::operator();\r
 222 \r
 223         double & operator() (size_t i, size_t j)\r
 224         {\r
 225                 return *gsl_matrix_ptr(m_matrix, i, j);\r
 226         }\r
 227 \r
 228         using base_type::get_gsl_matrix;\r
 229 \r
 230         gsl_matrix* get_gsl_matrix()\r
 231         {\r
 232                 return m_matrix;\r
 233         }\r
 234 \r
 235         VectorView row_view(size_t i)\r
 236         {\r
 237                 return VectorView(gsl_matrix_row(m_matrix, i));\r
 238         }\r
 239 \r
 240         VectorView column_view(size_t i)\r
 241         {\r
 242                 return VectorView(gsl_matrix_column(m_matrix, i));\r
 243         }\r
 244 \r
 245         void swap_rows(size_t i, size_t j)\r
 246         {\r
 247                  gsl_matrix_swap_rows(m_matrix, i, j);\r
 248         }\r
 249 \r
 250         void swap_columns(size_t i, size_t j)\r
 251         {\r
 252                 gsl_matrix_swap_columns(m_matrix, i, j);\r
 253         }\r
 254 \r
 255         MatrixImpl & transpose()\r
 256         {\r
 257             assert(columns() == rows());\r
 258                 gsl_matrix_transpose(m_matrix);\r
 259                 return (*this);\r
 260         }\r
 261 \r
 262         MatrixImpl & scale(double x)\r
 263         {\r
 264                 gsl_matrix_scale(m_matrix, x);\r
 265                 return (*this);\r
 266         }\r
 267 \r
 268         MatrixImpl & translate(double x)\r
 269         {\r
 270                  gsl_matrix_add_constant(m_matrix, x);\r
 271                  return (*this);\r
 272         }\r
 273 \r
 274         MatrixImpl & operator+=(base_type const& _matrix)\r
 275         {\r
 276                 gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());\r
 277                 return (*this);\r
 278         }\r
 279 \r
 280         MatrixImpl & operator-=(base_type const& _matrix)\r
 281         {\r
 282                 gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());\r
 283                 return (*this);\r
 284         }\r
 285 \r
 286 }; // end class MatrixImpl\r
 287 \r
 288 }  // end namespace detail\r
 289 \r
 290 \r
 291 using detail::operator==;\r
 292 using detail::operator<<;\r
 293 \r
 294 \r
 295 \r
 296 \r
 297 class Matrix: public detail::MatrixImpl\r
 298 {\r
 299   public:\r
 300         typedef detail::MatrixImpl base_type;\r
 301 \r
 302   public:\r
 303         // the matrix is not inizialized\r
 304         Matrix(size_t n1, size_t n2)\r
 305         {\r
 306                 m_rows = n1;\r
 307                 m_columns = n2;\r
 308                 m_matrix = gsl_matrix_alloc(n1, n2);\r
 309         }\r
 310 \r
 311         Matrix(size_t n1, size_t n2, double x)\r
 312         {\r
 313                 m_rows = n1;\r
 314                 m_columns = n2;\r
 315                 m_matrix = gsl_matrix_alloc(n1, n2);\r
 316                 gsl_matrix_set_all(m_matrix, x);\r
 317         }\r
 318 \r
 319         Matrix(Matrix const& _matrix)\r
 320         : base_type()\r
 321         {\r
 322                 m_rows = _matrix.rows();\r
 323                 m_columns = _matrix.columns();\r
 324                 m_matrix = gsl_matrix_alloc(rows(), columns());\r
 325                 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
 326         }\r
 327 \r
 328         explicit\r
 329         Matrix(base_type::base_type const& _matrix)\r
 330         {\r
 331                 m_rows = _matrix.rows();\r
 332                 m_columns = _matrix.columns();\r
 333                 m_matrix = gsl_matrix_alloc(rows(), columns());\r
 334                 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
 335         }\r
 336 \r
 337         Matrix & operator=(Matrix const& _matrix)\r
 338         {\r
 339                 assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );\r
 340                 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
 341                 return *this;\r
 342         }\r
 343 \r
 344         Matrix & operator=(base_type::base_type const& _matrix)\r
 345         {\r
 346                 assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );\r
 347                 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
 348                 return *this;\r
 349         }\r
 350 \r
 351         virtual ~Matrix()\r
 352         {\r
 353                 gsl_matrix_free(m_matrix);\r
 354         }\r
 355 \r
 356         Matrix & transpose()\r
 357         {\r
 358                 return static_cast<Matrix &>( base_type::transpose() );\r
 359         }\r
 360 \r
 361         Matrix & scale(double x)\r
 362         {\r
 363                 return static_cast<Matrix &>( base_type::scale(x) );\r
 364         }\r
 365 \r
 366         Matrix & translate(double x)\r
 367         {\r
 368                 return static_cast<Matrix &>( base_type::translate(x) );\r
 369         }\r
 370 \r
 371         Matrix & operator+=(base_type::base_type const& _matrix)\r
 372         {\r
 373                 return static_cast<Matrix &>( base_type::operator+=(_matrix) );\r
 374         }\r
 375 \r
 376         Matrix & operator-=(base_type::base_type const& _matrix)\r
 377         {\r
 378                 return static_cast<Matrix &>( base_type::operator-=(_matrix) );\r
 379         }\r
 380 \r
 381         friend\r
 382         void swap(Matrix & m1, Matrix & m2);\r
 383         friend\r
 384         void swap_any(Matrix & m1, Matrix & m2);\r
 385 \r
 386 };  // end class Matrix\r
 387 \r
 388 \r
 389 // warning! this operation invalidates any view of the passed matrix objects\r
 390 inline\r
 391 void swap(Matrix & m1, Matrix & m2)\r
 392 {\r
 393         assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );\r
 394         std::swap(m1.m_matrix, m2.m_matrix);\r
 395 }\r
 396 \r
 397 inline\r
 398 void swap_any(Matrix & m1, Matrix & m2)\r
 399 {\r
 400     std::swap(m1.m_matrix, m2.m_matrix);\r
 401     std::swap(m1.m_rows, m2.m_rows);\r
 402     std::swap(m1.m_columns, m2.m_columns);\r
 403 }\r
 404 \r
 405 \r
 406 \r
 407 class ConstMatrixView : public detail::BaseMatrixImpl\r
 408 {\r
 409   public:\r
 410         typedef detail::BaseMatrixImpl base_type;\r
 411 \r
 412   public:\r
 413         ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)\r
 414                 : m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )\r
 415         {\r
 416                 m_rows = n1;\r
 417                 m_columns = n2;\r
 418                 m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
 419         }\r
 420 \r
 421         ConstMatrixView(const ConstMatrixView & _matrix)\r
 422                 : base_type(),\r
 423                   m_matrix_view(_matrix.m_matrix_view)\r
 424         {\r
 425                 m_rows = _matrix.rows();\r
 426                 m_columns = _matrix.columns();\r
 427                 m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
 428         }\r
 429 \r
 430         ConstMatrixView(const base_type & _matrix)\r
 431                 : m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))\r
 432         {\r
 433                 m_rows = _matrix.rows();\r
 434                 m_columns = _matrix.columns();\r
 435                 m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
 436         }\r
 437 \r
 438   private:\r
 439         gsl_matrix_const_view m_matrix_view;\r
 440 \r
 441 };  // end class ConstMatrixView\r
 442 \r
 443 \r
 444 \r
 445 \r
 446 class MatrixView : public detail::MatrixImpl\r
 447 {\r
 448   public:\r
 449         typedef detail::MatrixImpl base_type;\r
 450 \r
 451   public:\r
 452         MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)\r
 453         {\r
 454                 m_rows = n1;\r
 455                 m_columns = n2;\r
 456                 m_matrix_view\r
 457                         = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);\r
 458                 m_matrix = &(m_matrix_view.matrix);\r
 459         }\r
 460 \r
 461         MatrixView(const MatrixView & _matrix)\r
 462         : base_type()\r
 463         {\r
 464                 m_rows = _matrix.rows();\r
 465                 m_columns = _matrix.columns();\r
 466                 m_matrix_view = _matrix.m_matrix_view;\r
 467                 m_matrix = &(m_matrix_view.matrix);\r
 468         }\r
 469 \r
 470         MatrixView(Matrix & _matrix)\r
 471         {\r
 472                 m_rows = _matrix.rows();\r
 473                 m_columns = _matrix.columns();\r
 474                 m_matrix_view\r
 475                         = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());\r
 476                 m_matrix = &(m_matrix_view.matrix);\r
 477         }\r
 478 \r
 479         MatrixView & operator=(MatrixView const& _matrix)\r
 480         {\r
 481                 assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );\r
 482                 gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);\r
 483                 return *this;\r
 484         }\r
 485 \r
 486         MatrixView & operator=(base_type::base_type const& _matrix)\r
 487         {\r
 488                 assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );\r
 489                 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
 490                 return *this;\r
 491         }\r
 492 \r
 493         MatrixView & transpose()\r
 494         {\r
 495                 return static_cast<MatrixView &>( base_type::transpose() );\r
 496         }\r
 497 \r
 498         MatrixView & scale(double x)\r
 499         {\r
 500                 return static_cast<MatrixView &>( base_type::scale(x) );\r
 501         }\r
 502 \r
 503         MatrixView & translate(double x)\r
 504         {\r
 505                 return static_cast<MatrixView &>( base_type::translate(x) );\r
 506         }\r
 507 \r
 508         MatrixView & operator+=(base_type::base_type const& _matrix)\r
 509         {\r
 510                 return static_cast<MatrixView &>( base_type::operator+=(_matrix) );\r
 511         }\r
 512 \r
 513         MatrixView & operator-=(base_type::base_type const& _matrix)\r
 514         {\r
 515                 return static_cast<MatrixView &>( base_type::operator-=(_matrix) );\r
 516         }\r
 517 \r
 518         friend\r
 519         void swap_view(MatrixView & m1, MatrixView & m2);\r
 520 \r
 521   private:\r
 522         gsl_matrix_view m_matrix_view;\r
 523 \r
 524 };  // end class MatrixView\r
 525 \r
 526 \r
 527 inline\r
 528 void swap_view(MatrixView & m1, MatrixView & m2)\r
 529 {\r
 530         assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );\r
 531         std::swap(m1.m_matrix_view, m2.m_matrix_view);\r
 532 }\r
 533 \r
 534 Vector operator*( detail::BaseMatrixImpl const& A,\r
 535                   detail::BaseVectorImpl const& v );\r
 536 \r
 537 Matrix operator*( detail::BaseMatrixImpl const& A,\r
 538                   detail::BaseMatrixImpl const& B );\r
 539 \r
 540 Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);\r
 541 \r
 542 } } // end namespaces\r
 543 \r
 544 #endif /*_NL_MATRIX_H_*/\r
 545 \r
 546 /*\r
 547   Local Variables:\r
 548   mode:c++\r
 549   c-file-style:"stroustrup"\r
 550   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
 551   indent-tabs-mode:nil\r
 552   fill-column:99\r
 553   End:\r
 554 */\r
 555 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r