src/2geom/numeric/vector.h

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