source: trunk/yat/utility/vector.h @ 783

Last change on this file since 783 was 783, checked in by Jari Häkkinen, 15 years ago

Fixes #155. Added asserts to check vectors in all non-const functions since no inline functions exists anymore.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.6 KB
Line 
1#ifndef _theplu_yat_utility_vector_
2#define _theplu_yat_utility_vector_
3
4// $Id: vector.h 783 2007-03-05 21:45:32Z jari $
5
6/*
7  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
8  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
9  Copyright (C) 2005 Jari Häkkinen, Peter Johansson, Markus Ringnér
10  Copyright (C) 2006, 2007 Jari Häkkinen
11
12  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
13
14  The yat library is free software; you can redistribute it and/or
15  modify it under the terms of the GNU General Public License as
16  published by the Free Software Foundation; either version 2 of the
17  License, or (at your option) any later version.
18
19  The yat library is distributed in the hope that it will be useful,
20  but WITHOUT ANY WARRANTY; without even the implied warranty of
21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22  General Public License for more details.
23
24  You should have received a copy of the GNU General Public License
25  along with this program; if not, write to the Free Software
26  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
27  02111-1307, USA.
28*/
29
30#include "Exception.h"
31
32#include <iostream>
33#include <vector>
34#include <utility>
35
36#include <gsl/gsl_vector.h>
37#include <gsl/gsl_sort_vector.h>
38
39namespace theplu {
40namespace yat {
41namespace utility {
42
43  class matrix;
44
45  /**
46     @brief This is the yat interface to GSL vector.
47
48     For time being 'double' is the only type supported.
49
50     \par File streams:
51     Reading and writing vectors to file streams are of course
52     supported. These are implemented without using GSL functionality,
53     and thus binary read and write to streams are not supported.
54
55     \par Vector views:
56     GSL vector views are supported and these are disguised as
57     ordinary utility::vectors. A support function is added,
58     utility::vector::isview(), that can be used to check if a vector
59     object is a view. Note that view vectors do not own the
60     underlying data, and a view is not valid if the vector owing the
61     data is deallocated.
62
63     \par
64     Currently there is no restriction on how a vector is used when
65     the vector is a const view into another vector or matrix. To
66     avoid unexpected runtime errors, the programmer must declare
67     const view vectors as 'const' in order to get compile time
68     diagnostics about improper use of a const view vector object. If
69     'const' is not used and the const view vector is used erroneously
70     (such as on the left hand side in assignments), the compiler will
71     not catch the error and runtime error will occur. assert(3) is
72     used to catch the runtime error during development. Example on
73     bad and proper use of const view vectors:
74     @code
75  const vector vm(13,1.0);
76  vector v1(vm,2,4);       // bad code! not const!
77  v1(0)=-123;              // accepted by compiler, runtime abort will occur
78                           // or catched by assert depending on compiler flags
79  const vector v2(vm,2,4); // proper code
80  v2(0)=-123;              // not acceptable for the compiler
81     @endcode
82
83     @note Missing support to underlying GSL vector features can in
84     principle be included to this class if requested by the user
85     community.
86  */
87
88  class vector
89  {
90  public:
91
92    /**
93       \brief The default constructor.
94    */
95    vector(void);
96
97    /**
98       \brief Allocates memory space for \a n elements, and sets all
99       elements to \a init_value.
100       
101       \throw GSL_error if memory allocation fails.
102    */
103    vector(size_t n, double init_value=0);
104
105    /**
106       \brief The copy constructor.
107
108       \note If the object to be copied is a vector view, the values
109       of the view will be copied, i.e. the view is not copied.
110
111       \throw A GSL_error is indirectly thrown if memory allocation
112       fails.
113    */
114    vector(const vector& other);
115
116    /**
117       \brief Vector view constructor.
118
119       Create a view of vector \a v, with starting index \a offset,
120       size \a n, and an optional \a stride.
121
122       A vector view can be used as any vector with the difference
123       that changes made to the view will also change the object that
124       is viewed. Also, using the copy constructor will create a new
125       vector object that is a copy of whatever is viewed. If a copy
126       of the view is needed then you should use this constructor to
127       obtain a copy.
128
129       \note If the object viewed by the view goes out of scope or is
130       deleted, the view becomes invalid and the result of further use
131       is undefined.
132
133       \throw GSL_error if a view cannot be set up.
134    */
135    vector(vector& v, size_t offset, size_t n, size_t stride=1);
136
137    /**
138       \brief Vector const view constructor.
139
140       Create a view of vector \a v, with starting index \a offset,
141       size \a n, and an optional \a stride.
142
143       A vector view can be used as any const vector. Using the copy
144       constructor will create a new vector object that is a copy of
145       whatever is viewed. If a copy of the view is needed then you
146       should use this constructor to obtain a copy.
147
148       \note If the object viewed by the view goes out of scope or is
149       deleted, the view becomes invalid and the result of further use
150       is undefined.
151
152       \throw GSL_error if a view cannot be set up.
153    */
154    vector(const vector& v, size_t offset, size_t n, size_t stride=1);
155
156    ///
157    /// Matrix row/column view constructor.
158    ///
159    /// Create a row/column vector view of matrix \a m, pointing at
160    /// row/column \a i. The parameter \a row is used to set whether
161    /// the view should be a row or column view. If \a row is set to
162    /// true, the view will be a row view (default behaviour), and,
163    /// naturally, a column view otherwise.
164    ///
165    /// A vector view can be used as any vector with the difference
166    /// that changes made to the view will also change the object that
167    /// is viewed. Also, using the copy constructor will create a new
168    /// vector object that is a copy of whatever is viewed. If a copy
169    /// of the view is needed then you should use the vector view
170    /// constructor to obtain a copy.
171    ///
172    /// @note If the object viewed by the view goes out of scope or is
173    /// deleted, the view becomes invalid and the result of further
174    /// use is undefined.
175    ///
176    vector(matrix& m, size_t i, bool row=true);
177
178    ///
179    /// Matrix row/column const view constructor.
180    ///
181    /// Create a row/column vector view of matrix \a m, pointing at
182    /// row/column \a i. The parameter \a row is used to set whether
183    /// the view should be a row or column view. If \a row is set to
184    /// true, the view will be a row view (default behaviour), and,
185    /// naturally, a column view otherwise.
186    ///
187    /// A const vector view can be used as any const vector. Using the
188    /// copy constructor will create a new vector object that is a
189    /// copy of whatever is viewed. If a copy of the view is needed
190    /// then you should use the vector view constructor to obtain a
191    /// copy.
192    ///
193    /// @note If the object viewed by the view goes out of scope or is
194    /// deleted, the view becomes invalid and the result of further
195    /// use is undefined.
196    ///
197    vector(const matrix& m, size_t i, bool row=true);
198
199    /**
200       \brief The istream constructor.
201
202       Either elements should be separated with white space characters
203       (default), or elements should be separated by the delimiter \a
204       sep. When delimiter \a sep is used empty elements are stored as
205       NaN's (except that empty lines are ignored). The end of input
206       to the vector is at end of file marker.
207
208       \throw GSL_error if memory allocation fails, IO_error if
209       unexpected input is found in the input stream.
210    */
211    explicit vector(std::istream &, char sep='\0')
212      throw(utility::IO_error, std::exception);
213
214    ///
215    /// The destructor.
216    ///
217    ~vector(void);
218
219    /**
220       \brief This function performs element-wise division, \f$ this_i =
221       this_i/other_i \; \forall i \f$.
222
223       \throw GSL_error if dimensions mis-match.
224    */
225    void div(const vector& other);
226
227    ///
228    /// @return A const pointer to the internal GSL vector,
229    ///
230    const gsl_vector* gsl_vector_p(void) const;
231
232    ///
233    /// @return A pointer to the internal GSL vector,
234    ///
235    gsl_vector* gsl_vector_p(void);
236
237    ///
238    /// Check if the vector object is a view (sub-vector) to another
239    /// vector.
240    ///
241    /// @return True if the object is a view, false othwerwise.
242    ///
243    bool isview(void) const;
244
245    /**
246       \brief This function performs element-wise multiplication, \f$
247       this_i = this_i * other_i \; \forall i \f$.
248
249       \throw GSL_error if dimensions mis-match.
250    */
251    void mul(const vector& other);
252
253    /**
254       \brief Reverse the order of elements in the vector.
255    */
256    void reverse(void);
257
258    /**
259       \brief Set element values to values in \a vec.
260
261       This function is needed for assignment of viewed elements.
262
263       \see const vector& operator=(const vector&)
264
265       \throw GSL_error if dimensions mis-match.
266    */
267    void set(const vector& vec);
268
269    ///
270    /// Set all elements to \a value.
271    ///
272    void set_all(const double& value);
273
274    ///
275    /// @return the number of elements in the vector.
276    ///
277    size_t size(void) const;
278
279    /**
280       \brief Exchange elements \a i and \a j.
281
282       \throw GSL_error if vector lengths differs.
283    */
284    void swap(size_t i, size_t j);
285
286    /**
287       \brief Element access operator.
288
289       \return Reference to element \a i.
290
291       \throw If GSL range checks are enabled in the underlying GSL
292       library a GSL_error exception is thrown if either index is out
293       of range.
294    */
295    double& operator()(size_t i);
296
297    /**
298       \brief Element access operator.
299
300       \return Const reference to element \a i.
301
302       \throw If GSL range checks are enabled in the underlying GSL
303       library a GSL_error exception is thrown if either index is out
304       of range.
305    */
306    const double& operator()(size_t i) const;
307
308    ///
309    /// Element access operator.
310    ///
311    /// @return Reference to element \a i.
312    ///
313    double& operator[](size_t i);
314
315    ///
316    /// Const element access operator.
317    ///
318    /// @return The value of element \a i.
319    ///
320    const double& operator[](size_t i) const;
321
322    ///
323    /// Comparison operator. Takes linear time.
324    ///
325    /// @return True if the sequence of the elements in the vectors
326    /// are element/wise equal.
327    ///
328    bool operator==(const vector&) const;
329
330    ///
331    /// @return The dot product.
332    ///
333    double operator*(const vector&) const;
334
335    /**
336       \brief The assignment operator.
337
338       There is no requirements on dimensions, i.e. the vector is
339       remapped in memory if necessary. This implies that in general
340       views cannot be assigned using this operator. Views will be
341       mutated into normal vectors. The only exception to this
342       behaviour on views is when self-assignemnt is done, since
343       self-assignment is ignored.
344
345       \return A const reference to the resulting vector.
346
347       \see void set(const vector&).
348
349       \throw A GSL_error is indirectly thrown if memory allocation
350       fails, or if dimensions mis-match.
351    */
352    const vector& operator=(const vector&);
353
354    /**
355       \brief Addition and assign operator. Vector addition, \f$
356       this_i = this_i + other_i \; \forall i \f$.
357
358       \return A const reference to the resulting vector.
359
360       \throw GSL_error if dimensions mis-match.
361    */
362    const vector& operator+=(const vector&);
363
364    /**
365       \brief Add a constant to a vector, \f$ this_i = this_i + d \;
366       \forall i \f$.
367
368       \return A const reference to the resulting vector.
369    */
370    const vector& operator+=(double d);
371
372    /**
373       \brief Subtract and assign operator. Vector subtraction, \f$
374       this_i = this_i - other_i \; \forall i \f$.
375
376       \return A const reference to the resulting vector.
377
378       \throw GSL_error if dimensions mis-match.
379    */
380    const vector& operator-=(const vector&);
381
382    /**
383       \brief Subtract a constant to a vector, \f$ this_i = this_i - d
384       \; \forall i \f$.
385
386       \return A const reference to the resulting vector.
387    */
388    const vector& operator-=(double d);
389
390    /**
391       \brief Multiply with scalar and assign operator, \f$ this_i =
392       this_i * d \; \forall i \f$.
393
394       \return A const reference to the resulting vector.
395    */
396    const vector& operator*=(double d);
397
398
399  private:
400
401    /**
402       \brief Create a new copy of the internal GSL vector.
403
404       Necessary memory for the new GSL vector is allocated and the
405       caller is responsible for freeing the allocated memory.
406
407       \return A pointer to a copy of the internal GSL vector.
408
409       \throw GSL_error if memory cannot be allocated for the new
410       copy, or if dimensions mis-match.
411    */
412    gsl_vector* create_gsl_vector_copy(void) const;
413
414    gsl_vector* v_;
415    const gsl_vector* v_const_;
416    gsl_vector_view* view_;
417    gsl_vector_const_view* view_const_;
418    // proxy_v_ is used to access the proper underlying v_ or v_const_
419    // in all const member functions. It is not used by design for
420    // non-const vector functions and operators. This is to make sure
421    // that runtime errors occur if a const vector is used in an
422    // inappropriate manner such as on left hand side in assignment
423    // (remember, v_ is null for const vector views).
424    const gsl_vector* proxy_v_;
425  };
426
427  /**
428     \brief Check if all elements of the vector are zero.
429
430     \return True if all elements in the vector is zero, false
431     othwerwise.
432  */
433  bool isnull(const vector&);
434
435  /**
436     \brief Get the maximum value of the vector.
437
438     \return The maximum value of the vector.
439  */
440  double max(const vector&);
441
442  /**
443     \brief Locate the maximum value in the vector.
444
445     \return The index to the maximum value of the vector.
446
447     \note Lower index has precedence.
448  */
449  size_t max_index(const vector&);
450
451  /**
452     \brief Get the minimum value of the vector.
453
454     \return The minimum value of the vector.
455  */
456  double min(const vector&);
457
458  /**
459     \brief Locate the minimum value in the vector.
460
461     \return The index to the minimum value of the vector.
462
463     \note Lower index has precedence.
464  */
465  size_t min_index(const vector&);
466
467  /**
468     \brief Get the minimum and maximim values of the vector.
469
470     \return The minimum and maximum values of the vector,
471     respectively.
472  */
473  std::pair<double,double> minmax(const vector&);
474
475  /**
476     \brief Locate the maximum and minumum element in the vector.
477
478     \return The indecies to the element with the minimum and maximum
479     values of the vector, respectively.
480
481     \note Lower index has precedence.
482  */
483  std::pair<size_t,size_t> minmax_index(const vector&);
484
485  /**
486     \brief Transforms a vector to a basis vector.
487
488     All elements are set to zero except the \a i-th element which is
489     set to one.
490  */
491  void set_basis(vector&, size_t i);
492
493  /**
494     Sort the elements in the vector.
495
496     \note Bug in GSL: if the vector contains one ore more NaNs the
497     call never returns.
498  */
499  void sort(vector&);
500
501  /**
502     \brief Calculate the sum of all vector elements.
503
504     \return The sum.
505  */
506  double sum(const vector&);
507
508  /**
509     \brief Swap vector elements by copying.
510
511     The two vectors must have the same length.
512
513     \throw GSL_error if vector lengths differs.
514  */
515  void swap(vector&, vector&);
516
517  /**
518     \brief The output operator for the vector class.
519  */
520  std::ostream& operator<<(std::ostream&, const vector&);
521
522}}} // of namespace utility, yat, and theplu
523
524#endif
Note: See TracBrowser for help on using the repository browser.