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

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

Addresses #193. vector now works as outlined here. Added some
functionality. Added a clone function that facilitates resizing of
vectors. clone is needed since assignement operator functionality is
changed.

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