src/2geom/numeric/vector.h

   1 /*\r
   2  * Vector, VectorView, ConstVectorView classes wrap the gsl vector 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_VECTOR_H_\r
  39 #define _NL_VECTOR_H_\r
  40 \r
  41 #include <cassert>\r
  42 #include <algorithm> // for std::swap\r
  43 #include <vector>\r
  44 #include <sstream>\r
  45 #include <string>\r
  46 \r
  47 \r
  48 #include <gsl/gsl_vector.h>\r
  49 #include <gsl/gsl_blas.h>\r
  50 \r
  51 \r
  52 namespace Geom { namespace NL {\r
  53 \r
  54 namespace detail\r
  55 {\r
  56 \r
  57 class BaseVectorImpl\r
  58 {\r
  59   public:\r
  60         double const& operator[](size_t i) const\r
  61         {\r
  62                 return *gsl_vector_const_ptr(m_vector, i);\r
  63         }\r
  64 \r
  65         const gsl_vector* get_gsl_vector() const\r
  66         {\r
  67                 return m_vector;\r
  68         }\r
  69         bool is_zero() const\r
  70         {\r
  71                 return gsl_vector_isnull(m_vector);\r
  72         }\r
  73 \r
  74         bool is_positive() const\r
  75         {\r
  76                 return gsl_vector_ispos(m_vector);\r
  77         }\r
  78 \r
  79         bool is_negative() const\r
  80         {\r
  81                 return gsl_vector_isneg(m_vector);\r
  82         }\r
  83 \r
  84         bool is_non_negative() const\r
  85         {\r
  86                 for ( size_t i = 0; i < size(); ++i )\r
  87                 {\r
  88                         if ( (*this)[i] < 0 ) return false;\r
  89                 }\r
  90                 return true;\r
  91         }\r
  92 \r
  93         double max() const\r
  94         {\r
  95                 return gsl_vector_max(m_vector);\r
  96         }\r
  97 \r
  98         double min() const\r
  99         {\r
 100                 return gsl_vector_min(m_vector);\r
 101         }\r
 102 \r
 103         size_t max_index() const\r
 104         {\r
 105                 return gsl_vector_max_index(m_vector);\r
 106         }\r
 107 \r
 108         size_t min_index() const\r
 109         {\r
 110                 return gsl_vector_min_index(m_vector);\r
 111         }\r
 112 \r
 113         size_t size() const\r
 114         {\r
 115                 return m_size;\r
 116         }\r
 117 \r
 118         std::string str() const;\r
 119 \r
 120         virtual ~BaseVectorImpl()\r
 121         {\r
 122         }\r
 123 \r
 124   protected:\r
 125         size_t m_size;\r
 126         gsl_vector* m_vector;\r
 127 \r
 128 };  // end class BaseVectorImpl\r
 129 \r
 130 \r
 131 inline\r
 132 bool operator== (BaseVectorImpl const& v1, BaseVectorImpl const& v2)\r
 133 {\r
 134         if (v1.size() != v2.size())     return false;\r
 135 \r
 136         for (size_t i = 0; i < v1.size(); ++i)\r
 137         {\r
 138                 if (v1[i] != v2[i])  return false;\r
 139         }\r
 140         return true;\r
 141 }\r
 142 \r
 143 template< class charT >\r
 144 inline\r
 145 std::basic_ostream<charT> &\r
 146 operator<< (std::basic_ostream<charT> & os, const BaseVectorImpl & _vector)\r
 147 {\r
 148         if (_vector.size() == 0 ) return os;\r
 149         os << "[" << _vector[0];\r
 150         for (unsigned int i = 1; i < _vector.size(); ++i)\r
 151         {\r
 152                 os << ", " << _vector[i];\r
 153         }\r
 154         os << "]";\r
 155         return os;\r
 156 }\r
 157 \r
 158 inline\r
 159 std::string BaseVectorImpl::str() const\r
 160 {\r
 161         std::ostringstream oss;\r
 162         oss << (*this);\r
 163         return oss.str();\r
 164 }\r
 165 \r
 166 inline\r
 167 double dot(BaseVectorImpl const& v1, BaseVectorImpl const& v2)\r
 168 {\r
 169     double result;\r
 170     gsl_blas_ddot(v1.get_gsl_vector(), v2.get_gsl_vector(), &result);\r
 171     return result;\r
 172 }\r
 173 \r
 174 \r
 175 class VectorImpl : public BaseVectorImpl\r
 176 {\r
 177   public:\r
 178         typedef BaseVectorImpl base_type;\r
 179 \r
 180   public:\r
 181         void set_all(double x)\r
 182         {\r
 183                 gsl_vector_set_all(m_vector, x);\r
 184         }\r
 185 \r
 186         void set_basis(size_t i)\r
 187         {\r
 188                 gsl_vector_set_basis(m_vector, i);\r
 189         }\r
 190 \r
 191         using base_type::operator[];\r
 192 \r
 193         double & operator[](size_t i)\r
 194         {\r
 195                 return *gsl_vector_ptr(m_vector, i);\r
 196         }\r
 197 \r
 198         using base_type::get_gsl_vector;\r
 199 \r
 200         gsl_vector* get_gsl_vector()\r
 201         {\r
 202                 return m_vector;\r
 203         }\r
 204 \r
 205         void swap_elements(size_t i, size_t j)\r
 206         {\r
 207                 gsl_vector_swap_elements(m_vector, i, j);\r
 208         }\r
 209 \r
 210         void reverse()\r
 211         {\r
 212                 gsl_vector_reverse(m_vector);\r
 213         }\r
 214 \r
 215         VectorImpl & scale(double x)\r
 216         {\r
 217                 gsl_vector_scale(m_vector, x);\r
 218                 return (*this);\r
 219         }\r
 220 \r
 221         VectorImpl & translate(double x)\r
 222         {\r
 223                 gsl_vector_add_constant(m_vector, x);\r
 224                 return (*this);\r
 225         }\r
 226 \r
 227         VectorImpl & operator+=(base_type const& _vector)\r
 228         {\r
 229                 gsl_vector_add(m_vector, _vector.get_gsl_vector());\r
 230                 return (*this);\r
 231         }\r
 232 \r
 233         VectorImpl & operator-=(base_type const& _vector)\r
 234         {\r
 235                 gsl_vector_sub(m_vector, _vector.get_gsl_vector());\r
 236                 return (*this);\r
 237         }\r
 238 \r
 239 };  // end class VectorImpl\r
 240 \r
 241 }  // end namespace detail\r
 242 \r
 243 \r
 244 using detail::operator==;\r
 245 using detail::operator<<;\r
 246 \r
 247 class Vector : public detail::VectorImpl\r
 248 {\r
 249   public:\r
 250         typedef detail::VectorImpl base_type;\r
 251 \r
 252   public:\r
 253         Vector(size_t n)\r
 254         {\r
 255                 m_size = n;\r
 256                 m_vector = gsl_vector_alloc(n);\r
 257         }\r
 258 \r
 259         Vector(size_t n, double x)\r
 260         {\r
 261                 m_size = n;\r
 262                 m_vector = gsl_vector_alloc(n);\r
 263                 gsl_vector_set_all(m_vector, x);\r
 264         }\r
 265 \r
 266         // create a vector with n elements all set to zero\r
 267         // but the i-th that is set to 1\r
 268         Vector(size_t n, size_t i)\r
 269         {\r
 270                 m_size = n;\r
 271                 m_vector = gsl_vector_alloc(n);\r
 272                 gsl_vector_set_basis(m_vector, i);\r
 273         }\r
 274 \r
 275         Vector(Vector const& _vector)\r
 276         : base_type()\r
 277         {\r
 278                 m_size = _vector.size();\r
 279                 m_vector = gsl_vector_alloc(size());\r
 280                 gsl_vector_memcpy(m_vector, _vector.m_vector);\r
 281         }\r
 282 \r
 283         explicit\r
 284         Vector(base_type::base_type const& _vector)\r
 285         {\r
 286                 m_size = _vector.size();\r
 287                 m_vector = gsl_vector_alloc(size());\r
 288                 gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
 289         }\r
 290 \r
 291         virtual ~Vector()\r
 292         {\r
 293                 gsl_vector_free(m_vector);\r
 294         }\r
 295 \r
 296 \r
 297         Vector & operator=(Vector const& _vector)\r
 298         {\r
 299                 assert( size() == _vector.size() );\r
 300                 gsl_vector_memcpy(m_vector, _vector.m_vector);\r
 301                 return (*this);\r
 302         }\r
 303 \r
 304         Vector & operator=(base_type::base_type const& _vector)\r
 305         {\r
 306                 assert( size() == _vector.size() );\r
 307                 gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
 308                 return (*this);\r
 309         }\r
 310 \r
 311         Vector & scale(double x)\r
 312         {\r
 313                 return static_cast<Vector&>( base_type::scale(x) );\r
 314         }\r
 315 \r
 316         Vector & translate(double x)\r
 317         {\r
 318                 return static_cast<Vector&>( base_type::translate(x) );\r
 319         }\r
 320 \r
 321         Vector & operator+=(base_type::base_type const& _vector)\r
 322         {\r
 323                 return static_cast<Vector&>( base_type::operator+=(_vector) );\r
 324         }\r
 325 \r
 326         Vector & operator-=(base_type::base_type const& _vector)\r
 327         {\r
 328                 return static_cast<Vector&>( base_type::operator-=(_vector) );\r
 329         }\r
 330 \r
 331         friend\r
 332         void swap(Vector & v1, Vector & v2);\r
 333         friend\r
 334         void swap_any(Vector & v1, Vector & v2);\r
 335 \r
 336 }; // end class Vector\r
 337 \r
 338 \r
 339 // warning! these operations invalidate any view of the passed vector objects\r
 340 inline\r
 341 void swap(Vector & v1, Vector & v2)\r
 342 {\r
 343         assert( v1.size() == v2.size() );\r
 344         std::swap(v1.m_vector, v2.m_vector);\r
 345 }\r
 346 \r
 347 inline\r
 348 void swap_any(Vector & v1, Vector & v2)\r
 349 {\r
 350     std::swap(v1.m_vector, v2.m_vector);\r
 351     std::swap(v1.m_size, v2.m_size);\r
 352 }\r
 353 \r
 354 \r
 355 class ConstVectorView : public detail::BaseVectorImpl\r
 356 {\r
 357   public:\r
 358         typedef detail::BaseVectorImpl base_type;\r
 359 \r
 360   public:\r
 361         ConstVectorView(const base_type & _vector, size_t n, size_t offset = 0)\r
 362                 : m_vector_view( gsl_vector_const_subvector(_vector.get_gsl_vector(), offset, n) )\r
 363         {\r
 364                 m_size = n;\r
 365                 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
 366         }\r
 367 \r
 368         ConstVectorView(const base_type & _vector, size_t n, size_t offset , size_t stride)\r
 369                 : m_vector_view( gsl_vector_const_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n) )\r
 370         {\r
 371                 m_size = n;\r
 372                 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
 373         }\r
 374 \r
 375     ConstVectorView(const double* _vector, size_t n, size_t offset = 0)\r
 376         : m_vector_view( gsl_vector_const_view_array(_vector + offset, n) )\r
 377     {\r
 378         m_size = n;\r
 379         m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
 380     }\r
 381 \r
 382     ConstVectorView(const double* _vector, size_t n, size_t offset, size_t stride)\r
 383         : m_vector_view( gsl_vector_const_view_array_with_stride(_vector + offset, stride, n) )\r
 384     {\r
 385         m_size = n;\r
 386         m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
 387     }\r
 388 \r
 389         explicit\r
 390         ConstVectorView(gsl_vector_const_view  _gsl_vector_view)\r
 391                 : m_vector_view(_gsl_vector_view)\r
 392         {\r
 393                 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
 394                 m_size = m_vector->size;\r
 395         }\r
 396 \r
 397     explicit\r
 398     ConstVectorView(const std::vector<double>&  _vector)\r
 399         : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) )\r
 400     {\r
 401         m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
 402         m_size = _vector.size();\r
 403     }\r
 404 \r
 405         ConstVectorView(const ConstVectorView & _vector)\r
 406                 : base_type(),\r
 407                   m_vector_view(_vector.m_vector_view)\r
 408         {\r
 409                 m_size = _vector.size();\r
 410                 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
 411         }\r
 412 \r
 413         ConstVectorView(const base_type & _vector)\r
 414                 : m_vector_view(gsl_vector_const_subvector(_vector.get_gsl_vector(), 0, _vector.size()))\r
 415         {\r
 416                 m_size = _vector.size();\r
 417                 m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
 418         }\r
 419 \r
 420   private:\r
 421         gsl_vector_const_view m_vector_view;\r
 422 \r
 423 }; // end class  ConstVectorView\r
 424 \r
 425 \r
 426 \r
 427 \r
 428 class VectorView : public detail::VectorImpl\r
 429 {\r
 430   public:\r
 431         typedef detail::VectorImpl base_type;\r
 432 \r
 433   public:\r
 434         VectorView(base_type & _vector, size_t n, size_t offset = 0, size_t stride = 1)\r
 435         {\r
 436                 m_size = n;\r
 437                 if (stride == 1)\r
 438                 {\r
 439                         m_vector_view\r
 440                                 = gsl_vector_subvector(_vector.get_gsl_vector(), offset, n);\r
 441                         m_vector = &(m_vector_view.vector);\r
 442                 }\r
 443                 else\r
 444                 {\r
 445                         m_vector_view\r
 446                                 = gsl_vector_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n);\r
 447                         m_vector = &(m_vector_view.vector);\r
 448                 }\r
 449         }\r
 450 \r
 451     VectorView(double* _vector, size_t n, size_t offset = 0, size_t stride = 1)\r
 452     {\r
 453         m_size = n;\r
 454         if (stride == 1)\r
 455         {\r
 456             m_vector_view\r
 457                 = gsl_vector_view_array(_vector + offset, n);\r
 458             m_vector = &(m_vector_view.vector);\r
 459         }\r
 460         else\r
 461         {\r
 462             m_vector_view\r
 463                 = gsl_vector_view_array_with_stride(_vector + offset, stride, n);\r
 464             m_vector = &(m_vector_view.vector);\r
 465         }\r
 466 \r
 467     }\r
 468 \r
 469         VectorView(const VectorView & _vector)\r
 470         : base_type()\r
 471         {\r
 472                 m_size = _vector.size();\r
 473                 m_vector_view = _vector.m_vector_view;\r
 474                 m_vector = &(m_vector_view.vector);\r
 475         }\r
 476 \r
 477         VectorView(Vector & _vector)\r
 478         {\r
 479                 m_size = _vector.size();\r
 480                 m_vector_view = gsl_vector_subvector(_vector.get_gsl_vector(), 0, size());\r
 481                 m_vector = &(m_vector_view.vector);\r
 482         }\r
 483 \r
 484         explicit\r
 485         VectorView(gsl_vector_view _gsl_vector_view)\r
 486                 : m_vector_view(_gsl_vector_view)\r
 487         {\r
 488                 m_vector = &(m_vector_view.vector);\r
 489                 m_size = m_vector->size;\r
 490         }\r
 491 \r
 492         explicit\r
 493         VectorView(std::vector<double> & _vector)\r
 494         {\r
 495             m_size = _vector.size();\r
 496             m_vector_view = gsl_vector_view_array(&(_vector[0]), _vector.size());\r
 497             m_vector = &(m_vector_view.vector);\r
 498         }\r
 499 \r
 500         VectorView & operator=(VectorView const& _vector)\r
 501         {\r
 502                 assert( size() == _vector.size() );\r
 503                 gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
 504                 return (*this);\r
 505         }\r
 506 \r
 507         VectorView & operator=(base_type::base_type const& _vector)\r
 508         {\r
 509                 assert( size() == _vector.size() );\r
 510                 gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
 511                 return (*this);\r
 512         }\r
 513 \r
 514         VectorView & scale(double x)\r
 515         {\r
 516                 return static_cast<VectorView&>( base_type::scale(x) );\r
 517         }\r
 518 \r
 519         VectorView & translate(double x)\r
 520         {\r
 521                 return static_cast<VectorView&>( base_type::translate(x) );\r
 522         }\r
 523 \r
 524         VectorView & operator+=(base_type::base_type const& _vector)\r
 525         {\r
 526                 return static_cast<VectorView&>( base_type::operator+=(_vector) );\r
 527         }\r
 528 \r
 529         VectorView & operator-=(base_type::base_type const& _vector)\r
 530         {\r
 531                 return static_cast<VectorView&>( base_type::operator-=(_vector) );\r
 532         }\r
 533 \r
 534         friend\r
 535         void swap_view(VectorView & v1, VectorView & v2);\r
 536 \r
 537   private:\r
 538         gsl_vector_view m_vector_view;\r
 539 \r
 540 }; // end class VectorView\r
 541 \r
 542 \r
 543 inline\r
 544 void swap_view(VectorView & v1, VectorView & v2)\r
 545 {\r
 546         assert( v1.size() == v2.size() );\r
 547         std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too\r
 548 }\r
 549 \r
 550 \r
 551 } } // end namespaces\r
 552 \r
 553 \r
 554 #endif /*_NL_VECTOR_H_*/\r
 555 \r
 556 /*\r
 557   Local Variables:\r
 558   mode:c++\r
 559   c-file-style:"stroustrup"\r
 560   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
 561   indent-tabs-mode:nil\r
 562   fill-column:99\r
 563   End:\r
 564 */\r
 565 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r