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

Last change on this file since 825 was 825, checked in by Peter, 15 years ago

removed minmax functions and moved shuffle into vector.*

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