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