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

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

Fixes #195.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.2 KB
Line 
1#ifndef _theplu_yat_utility_vector_
2#define _theplu_yat_utility_vector_
3
4// $Id: vector.h 782 2007-03-05 21:17: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, 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
66     matrix. To avoid unexpected runtime errors, the programmer must
67     declare const view vectors as const in order to get compile time
68     diagnostics about improper use of a const view vector object. An
69     example of improper use of a const view vector is assignment of
70     values to an element in the object. If the const view object was
71     declared const the assigment will be caught by the compiler
72     whereas it will cause a (un-catchable) runtime error. Example:
73     @code
74  const vector vm(13,1.0);
75  vector v1(vm,2,4);       // bad code! not const!
76  v1(0)=-123;              // accepted by compiler, runtime abort will occur
77  const vector v2(vm,2,4); // proper code
78  v2(0)=-123;              // not acceptable for the compiler
79     @endcode
80
81     @note Missing support to underlying GSL vector features can in
82     principle be included to this class if requested by the user
83     community.
84  */
85
86  class vector
87  {
88  public:
89
90    /**
91       \brief The default constructor.
92    */
93    vector(void);
94
95    /**
96       \brief Allocates memory space for \a n elements, and sets all
97       elements to \a init_value.
98       
99       \throw GSL_error if memory allocation fails.
100    */
101    vector(size_t n, double init_value=0);
102
103    /**
104       \brief The copy constructor.
105
106       \note If the object to be copied is a vector view, the values
107       of the view will be copied, i.e. the view is not copied.
108
109       \throw A GSL_error is indirectly thrown if memory allocation
110       fails.
111    */
112    vector(const vector& other);
113
114    /**
115       \brief Vector view constructor.
116
117       Create a view of vector \a v, with starting index \a offset,
118       size \a n, and an optional \a stride.
119
120       A vector view can be used as any vector with the difference
121       that changes made to the view will also change the object that
122       is viewed. Also, using the copy constructor will create a new
123       vector object that is a copy of whatever is viewed. If a copy
124       of the view is needed then you should use this constructor to
125       obtain a copy.
126
127       \note If the object viewed by the view goes out of scope or is
128       deleted, the view becomes invalid and the result of further use
129       is undefined.
130
131       \throw GSL_error if a view cannot be set up.
132    */
133    vector(vector& v, size_t offset, size_t n, size_t stride=1);
134
135    /**
136       \brief Vector const view constructor.
137
138       Create a view of vector \a v, with starting index \a offset,
139       size \a n, and an optional \a stride.
140
141       A vector view can be used as any const vector. Using the copy
142       constructor will create a new vector object that is a copy of
143       whatever is viewed. If a copy of the view is needed then you
144       should use this constructor to obtain a copy.
145
146       \note If the object viewed by the view goes out of scope or is
147       deleted, the view becomes invalid and the result of further use
148       is undefined.
149
150       \throw GSL_error if a view cannot be set up.
151    */
152    vector(const vector& v, size_t offset, size_t n, size_t stride=1);
153
154    ///
155    /// Matrix row/column view constructor.
156    ///
157    /// Create a row/column vector view of matrix \a m, pointing at
158    /// row/column \a i. The parameter \a row is used to set whether
159    /// the view should be a row or column view. If \a row is set to
160    /// true, the view will be a row view (default behaviour), and,
161    /// naturally, a column view otherwise.
162    ///
163    /// A vector view can be used as any vector with the difference
164    /// that changes made to the view will also change the object that
165    /// is viewed. Also, using the copy constructor will create a new
166    /// vector object that is a copy of whatever is viewed. If a copy
167    /// of the view is needed then you should use the vector view
168    /// constructor to obtain a copy.
169    ///
170    /// @note If the object viewed by the view goes out of scope or is
171    /// deleted, the view becomes invalid and the result of further
172    /// use is undefined.
173    ///
174    vector(matrix& m, size_t i, bool row=true);
175
176    ///
177    /// Matrix row/column const view constructor.
178    ///
179    /// Create a row/column vector view of matrix \a m, pointing at
180    /// row/column \a i. The parameter \a row is used to set whether
181    /// the view should be a row or column view. If \a row is set to
182    /// true, the view will be a row view (default behaviour), and,
183    /// naturally, a column view otherwise.
184    ///
185    /// A const vector view can be used as any const vector. Using the
186    /// copy constructor will create a new vector object that is a
187    /// copy of whatever is viewed. If a copy of the view is needed
188    /// then you should use the vector view constructor to obtain a
189    /// copy.
190    ///
191    /// @note If the object viewed by the view goes out of scope or is
192    /// deleted, the view becomes invalid and the result of further
193    /// use is undefined.
194    ///
195    vector(const matrix& m, size_t i, bool row=true);
196
197    /**
198       \brief The istream constructor.
199
200       Either elements should be separated with white space characters
201       (default), or elements should be separated by the delimiter \a
202       sep. When delimiter \a sep is used empty elements are stored as
203       NaN's (except that empty lines are ignored). The end of input
204       to the vector is at end of file marker.
205
206       \throw GSL_error if memory allocation fails, IO_error if
207       unexpected input is found in the input stream.
208    */
209    explicit vector(std::istream &, char sep='\0')
210      throw(utility::IO_error, std::exception);
211
212    ///
213    /// The destructor.
214    ///
215    ~vector(void);
216
217    /**
218       \brief This function performs element-wise division, \f$ this_i =
219       this_i/other_i \; \forall i \f$.
220
221       \throw GSL_error if dimensions mis-match.
222    */
223    void div(const vector& other);
224
225    ///
226    /// @return A const pointer to the internal GSL vector,
227    ///
228    const gsl_vector* gsl_vector_p(void) const;
229
230    ///
231    /// @return A pointer to the internal GSL vector,
232    ///
233    gsl_vector* gsl_vector_p(void);
234
235    ///
236    /// Check if the vector object is a view (sub-vector) to another
237    /// vector.
238    ///
239    /// @return True if the object is a view, false othwerwise.
240    ///
241    bool isview(void) const;
242
243    /**
244       \brief This function performs element-wise multiplication, \f$
245       this_i = this_i * other_i \; \forall i \f$.
246
247       \throw GSL_error if dimensions mis-match.
248    */
249    void mul(const vector& other);
250
251    /**
252       \brief Reverse the order of elements in the vector.
253    */
254    void reverse(void);
255
256    /**
257       \brief Set element values to values in \a vec.
258
259       This function is needed for assignment of viewed elements.
260
261       \see const vector& operator=(const vector&)
262
263       \throw GSL_error if dimensions mis-match.
264    */
265    void set(const vector& vec);
266
267    ///
268    /// Set all elements to \a value.
269    ///
270    void set_all(const double& value);
271
272    ///
273    /// @return the number of elements in the vector.
274    ///
275    size_t size(void) const;
276
277    /**
278       \brief Exchange elements \a i and \a j.
279
280       \throw GSL_error if vector lengths differs.
281    */
282    void swap(size_t i, size_t j);
283
284    /**
285       \brief Element access operator.
286
287       \return Reference to element \a i.
288
289       \throw If GSL range checks are enabled in the underlying GSL
290       library a GSL_error exception is thrown if either index is out
291       of range.
292    */
293    double& operator()(size_t i);
294
295    /**
296       \brief Element access operator.
297
298       \return Const reference to element \a i.
299
300       \throw If GSL range checks are enabled in the underlying GSL
301       library a GSL_error exception is thrown if either index is out
302       of range.
303    */
304    const double& operator()(size_t i) const;
305
306    ///
307    /// Element access operator.
308    ///
309    /// @return Reference to element \a i.
310    ///
311    double& operator[](size_t i);
312
313    ///
314    /// Const element access operator.
315    ///
316    /// @return The value of element \a i.
317    ///
318    const double& operator[](size_t i) const;
319
320    ///
321    /// Comparison operator. Takes linear time.
322    ///
323    /// @return True if the sequence of the elements in the vectors
324    /// are element/wise equal.
325    ///
326    bool operator==(const vector&) const;
327
328    ///
329    /// @return The dot product.
330    ///
331    double operator*(const vector&) const;
332
333    /**
334       \brief The assignment operator.
335
336       There is no requirements on dimensions, i.e. the vector is
337       remapped in memory if necessary. This implies that in general
338       views cannot be assigned using this operator. Views will be
339       mutated into normal vectors. The only exception to this
340       behaviour on views is when self-assignemnt is done, since
341       self-assignment is ignored.
342
343       \return A const reference to the resulting vector.
344
345       \see void set(const vector&).
346
347       \throw A GSL_error is indirectly thrown if memory allocation
348       fails, or if dimensions mis-match.
349    */
350    const vector& operator=(const vector&);
351
352    /**
353       \brief Addition and assign operator. Vector addition, \f$
354       this_i = this_i + other_i \; \forall i \f$.
355
356       \return A const reference to the resulting vector.
357
358       \throw GSL_error if dimensions mis-match.
359    */
360    const vector& operator+=(const vector&);
361
362    /**
363       \brief Add a constant to a vector, \f$ this_i = this_i + d \;
364       \forall i \f$.
365
366       \return A const reference to the resulting vector.
367    */
368    const vector& operator+=(double d);
369
370    /**
371       \brief Subtract and assign operator. Vector subtraction, \f$
372       this_i = this_i - other_i \; \forall i \f$.
373
374       \return A const reference to the resulting vector.
375
376       \throw GSL_error if dimensions mis-match.
377    */
378    const vector& operator-=(const vector&);
379
380    /**
381       \brief Subtract a constant to a vector, \f$ this_i = this_i - d
382       \; \forall i \f$.
383
384       \return A const reference to the resulting vector.
385    */
386    const vector& operator-=(double d);
387
388    /**
389       \brief Multiply with scalar and assign operator, \f$ this_i =
390       this_i * d \; \forall i \f$.
391
392       \return A const reference to the resulting vector.
393    */
394    const vector& operator*=(double d);
395
396
397  private:
398
399    /**
400       \brief Create a new copy of the internal GSL vector.
401
402       Necessary memory for the new GSL vector is allocated and the
403       caller is responsible for freeing the allocated memory.
404
405       \return A pointer to a copy of the internal GSL vector.
406
407       \throw GSL_error if memory cannot be allocated for the new
408       copy, or if dimensions mis-match.
409    */
410    gsl_vector* create_gsl_vector_copy(void) const;
411
412    gsl_vector* v_;
413    const gsl_vector* v_const_;
414    gsl_vector_view* view_;
415    gsl_vector_const_view* view_const_;
416    // proxy_v_ is used to access the proper underlying v_ or v_const_
417    // in all const member functions.
418    const gsl_vector* proxy_v_;
419  };
420
421  /**
422     \brief Check if all elements of the vector are zero.
423
424     \return True if all elements in the vector is zero, false
425     othwerwise.
426  */
427  bool isnull(const vector&);
428
429  /**
430     \brief Get the maximum value of the vector.
431
432     \return The maximum value of the vector.
433  */
434  double max(const vector&);
435
436  /**
437     \brief Locate the maximum value in the vector.
438
439     \return The index to the maximum value of the vector.
440
441     \note Lower index has precedence.
442  */
443  size_t max_index(const vector&);
444
445  /**
446     \brief Get the minimum value of the vector.
447
448     \return The minimum value of the vector.
449  */
450  double min(const vector&);
451
452  /**
453     \brief Locate the minimum value in the vector.
454
455     \return The index to the minimum value of the vector.
456
457     \note Lower index has precedence.
458  */
459  size_t min_index(const vector&);
460
461  /**
462     \brief Get the minimum and maximim values of the vector.
463
464     \return The minimum and maximum values of the vector,
465     respectively.
466  */
467  std::pair<double,double> minmax(const vector&);
468
469  /**
470     \brief Locate the maximum and minumum element in the vector.
471
472     \return The indecies to the element with the minimum and maximum
473     values of the vector, respectively.
474
475     \note Lower index has precedence.
476  */
477  std::pair<size_t,size_t> minmax_index(const vector&);
478
479  /**
480     \brief Transforms a vector to a basis vector.
481
482     All elements are set to zero except the \a i-th element which is
483     set to one.
484  */
485  void set_basis(vector&, size_t i);
486
487  /**
488     Sort the elements in the vector.
489
490     \note Bug in GSL: if the vector contains one ore more NaNs the
491     call never returns.
492  */
493  void sort(vector&);
494
495  /**
496     \brief Calculate the sum of all vector elements.
497
498     \return The sum.
499  */
500  double sum(const vector&);
501
502  /**
503     \brief Swap vector elements by copying.
504
505     The two vectors must have the same length.
506
507     \throw GSL_error if vector lengths differs.
508  */
509  void swap(vector&, vector&);
510
511  /**
512     \brief The output operator for the vector class.
513  */
514  std::ostream& operator<<(std::ostream&, const vector&);
515
516}}} // of namespace utility, yat, and theplu
517
518#endif
Note: See TracBrowser for help on using the repository browser.