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

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

Addresses #193. Fixing doxygen docs, and operator=.

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