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

Last change on this file since 865 was 865, checked in by Peter, 14 years ago

changing URL to http://trac.thep.lu.se/trac/yat

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