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

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

Fixed spelling.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.0 KB
Line 
1#ifndef _theplu_yat_utility_vector_
2#define _theplu_yat_utility_vector_
3
4// $Id: vector.h 809 2007-03-15 22:44:36Z 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 Resize vector
271    ///
272    /// All elements are set to @a init_value.
273    ///
274    void resize(size_t, double init_value=0);
275
276    /**
277       \brief Reverse the order of elements in the vector.
278    */
279    void reverse(void);
280
281    /**
282       \brief Set element values to values in \a vec.
283
284       This function is needed for assignment of viewed elements.
285
286       \see const vector& operator=(const vector&)
287
288       \throw GSL_error if dimensions mis-match.
289    */
290    void set(const vector& vec);
291
292    ///
293    /// Set all elements to \a value.
294    ///
295    void set_all(const double& value);
296
297    ///
298    /// @return the number of elements in the vector.
299    ///
300    size_t size(void) const;
301
302    /**
303       \brief Exchange elements \a i and \a j.
304
305       \throw GSL_error if vector lengths differs.
306    */
307    void swap(size_t i, size_t j);
308
309    /**
310       \brief Element access operator.
311
312       \return Reference to element \a i.
313
314       \throw If GSL range checks are enabled in the underlying GSL
315       library a GSL_error exception is thrown if either index is out
316       of range.
317    */
318    double& operator()(size_t i);
319
320    /**
321       \brief Element access operator.
322
323       \return Const reference to element \a i.
324
325       \throw If GSL range checks are enabled in the underlying GSL
326       library a GSL_error exception is thrown if either index is out
327       of range.
328    */
329    const double& operator()(size_t i) const;
330
331    ///
332    /// Element access operator.
333    ///
334    /// @return Reference to element \a i.
335    ///
336    double& operator[](size_t i);
337
338    ///
339    /// Const element access operator.
340    ///
341    /// @return The value of element \a i.
342    ///
343    const double& operator[](size_t i) const;
344
345    /**
346       \brief Comparison operator. Takes linear time.
347
348       Checks are performed with exact matching, i.e., rounding off
349       effects may destroy comparison. Use the equal function for
350       comparing elements within a user defined precision.
351
352       \return True if all elements are equal otherwise false.
353
354       \see equal(const vector&, const double precision=0)
355    */
356    bool operator==(const vector&) const;
357
358    /**
359       \brief Comparison operator. Takes linear time.
360
361       Checks are performed with exact matching, i.e., rounding off
362       effects may destroy comparison. Use the equal function for
363       comparing elements within a user defined precision.
364
365       \return False if all elements are equal otherwise true.
366
367       \see equal(const vector&, const double precision=0)
368    */
369    bool operator!=(const vector&) const;
370
371    ///
372    /// @return The dot product.
373    ///
374    double operator*(const vector&) const;
375
376    /**
377       \brief The assignment operator.
378
379       Dimensions of the vectors must match. If the LHS vector is a
380       view, the underlying data will be changed.
381
382       \return A const reference to the resulting vector.
383
384       \see void set(const vector&).
385
386       \throw GSL_error if dimensions mis-match.
387    */
388    const vector& operator=(const vector&);
389
390    /**
391       \brief Addition and assign operator. Vector addition, \f$
392       this_i = this_i + other_i \; \forall i \f$.
393
394       \return A const reference to the resulting vector.
395
396       \throw GSL_error if dimensions mis-match.
397    */
398    const vector& operator+=(const vector&);
399
400    /**
401       \brief Add a constant to a vector, \f$ this_i = this_i + d \;
402       \forall i \f$.
403
404       \return A const reference to the resulting vector.
405    */
406    const vector& operator+=(double d);
407
408    /**
409       \brief Subtract and assign operator. Vector subtraction, \f$
410       this_i = this_i - other_i \; \forall i \f$.
411
412       \return A const reference to the resulting vector.
413
414       \throw GSL_error if dimensions mis-match.
415    */
416    const vector& operator-=(const vector&);
417
418    /**
419       \brief Subtract a constant to a vector, \f$ this_i = this_i - d
420       \; \forall i \f$.
421
422       \return A const reference to the resulting vector.
423    */
424    const vector& operator-=(double d);
425
426    /**
427       \brief Multiply with scalar and assign operator, \f$ this_i =
428       this_i * d \; \forall i \f$.
429
430       \return A const reference to the resulting vector.
431    */
432    const vector& operator*=(double d);
433
434
435  private:
436
437    /**
438       \brief Create a new copy of the internal GSL vector.
439
440       Necessary memory for the new GSL vector is allocated and the
441       caller is responsible for freeing the allocated memory.
442
443       \return A pointer to a copy of the internal GSL vector.
444
445       \throw GSL_error if memory cannot be allocated for the new
446       copy, or if dimensions mis-match.
447    */
448    gsl_vector* create_gsl_vector_copy(void) const;
449
450    gsl_vector* v_;
451    gsl_vector_view* view_;
452    const gsl_vector_const_view* view_const_;
453    // proxy_v_ is used to access the proper underlying gsl_vector
454    // in all const member functions. It is not used by design for
455    // non-const vector functions and operators. This is to make sure
456    // that runtime errors occur if a const vector is used in an
457    // inappropriate manner such as on left hand side in assignment
458    // (remember, v_ is null for const vector views).
459    const gsl_vector* proxy_v_;
460  };
461
462  /**
463     \brief Check if all elements of the vector are zero.
464
465     \return True if all elements in the vector is zero, false
466     othwerwise.
467  */
468  bool isnull(const vector&);
469
470  /**
471     \brief Get the maximum value of the vector.
472
473     \return The maximum value of the vector.
474  */
475  double max(const vector&);
476
477  /**
478     \brief Locate the maximum value in the vector.
479
480     \return The index to the maximum value of the vector.
481
482     \note Lower index has precedence.
483  */
484  size_t max_index(const vector&);
485
486  /**
487     \brief Get the minimum value of the vector.
488
489     \return The minimum value of the vector.
490  */
491  double min(const vector&);
492
493  /**
494     \brief Locate the minimum value in the vector.
495
496     \return The index to the minimum value of the vector.
497
498     \note Lower index has precedence.
499  */
500  size_t min_index(const vector&);
501
502  /**
503     \brief Get the minimum and maximim values of the vector.
504
505     \return The minimum and maximum values of the vector,
506     respectively.
507  */
508  std::pair<double,double> minmax(const vector&);
509
510  /**
511     \brief Locate the maximum and minumum element in the vector.
512
513     \return The indecies to the element with the minimum and maximum
514     values of the vector, respectively.
515
516     \note Lower index has precedence.
517  */
518  std::pair<size_t,size_t> minmax_index(const vector&);
519
520  /**
521     \brief Create a vector \a flag indicating NaN's in another vector
522     \a templat.
523
524     The \a flag vector is changed to contain 1's and 0's only. A 1
525     means that the corresponding element in the \a templat vector is
526     valid and a zero means that the corresponding element is a NaN.
527
528     \note Space for vector \a flag is reallocated to fit the size of
529     vector \a templat if sizes mismatch.
530
531     \return True if the \a templat vector contains at least one NaN.
532  */
533  bool nan(const matrix& templat, matrix& flag);
534
535  /**
536     \brief Transforms a vector to a basis vector.
537
538     All elements are set to zero except the \a i-th element which is
539     set to one.
540  */
541  void set_basis(vector&, size_t i);
542
543  /**
544     Sort the elements in the vector.
545
546     \note Bug in GSL: if the vector contains one ore more NaNs the
547     call never returns.
548  */
549  void sort(vector&);
550
551  /**
552     \brief Calculate the sum of all vector elements.
553
554     \return The sum.
555  */
556  double sum(const vector&);
557
558  /**
559     \brief Swap vector elements by copying.
560
561     The two vectors must have the same length.
562
563     \throw GSL_error if vector lengths differs.
564  */
565  void swap(vector&, vector&);
566
567  /**
568     \brief The output operator for the vector class.
569  */
570  std::ostream& operator<<(std::ostream&, const vector&);
571
572}}} // of namespace utility, yat, and theplu
573
574#endif
Note: See TracBrowser for help on using the repository browser.