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

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

merged constIterator into Iterator and added an extra template parameter. There are problems with conversion from iterator to const_iterator

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.9 KB
Line 
1#ifndef _theplu_yat_utility_vector_
2#define _theplu_yat_utility_vector_
3
4// $Id: vector.h 881 2007-09-22 20:22:36Z 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 "constIterator.h"
32#include "Exception.h"
33#include "Iterator.h"
34
35#include <iostream>
36#include <vector>
37#include <utility>
38
39#include <gsl/gsl_vector.h>
40#include <gsl/gsl_sort_vector.h>
41
42namespace theplu {
43namespace yat {
44namespace utility {
45
46  class matrix;
47
48  /**
49     @brief This is the yat interface to GSL vector.
50
51     For time being 'double' is the only type supported.
52
53     \par File streams:
54     Reading and writing vectors to file streams are of course
55     supported. These are implemented without using GSL functionality,
56     and thus binary read and write to streams are not supported.
57
58     \par Vector views:
59     GSL vector views are supported and these are disguised as
60     ordinary utility::vectors. A support function is added,
61     utility::vector::isview(), that can be used to check if a vector
62     object is a view. Note that view vectors do not own the
63     underlying data, and a view is not valid if the vector/matrix
64     owing the data is deallocated.
65
66     \par
67     Currently there is no restriction on how a vector is used when
68     the vector is a const view into another vector or matrix. To
69     avoid unexpected runtime errors, the programmer must declare
70     const view vectors as 'const' in order to get compile time
71     diagnostics about improper use of a const view vector object. If
72     'const' is not used and the const view vector is used erroneously
73     (such as on the left hand side in assignments), the compiler will
74     not catch the error and runtime error will occur. assert(3) is
75     used to catch the runtime error during development. Example on
76     bad and proper use of const view vectors:
77     @code
78  const vector vm(13,1.0);
79  vector v1(vm,2,4);       // bad code! not const!
80  v1(0)=-123;              // accepted by compiler, runtime abort will occur
81                           // or catched by assert depending on compiler flags
82  const vector v2(vm,2,4); // proper code
83  v2(0)=-123;              // not acceptable for the compiler
84     @endcode
85  */
86
87  class vector
88  {
89  public:
90
91    typedef Iterator<double&, vector> iterator;
92    typedef Iterator<const double, const vector> const_iterator;
93
94    /**
95       \brief The default constructor.
96    */
97    vector(void);
98
99    /**
100       \brief Allocates memory space for \a n elements, and sets all
101       elements to \a init_value.
102       
103       \throw GSL_error if memory allocation fails.
104    */
105    vector(size_t n, double init_value=0);
106
107    /**
108       \brief The copy constructor.
109
110       \note If the object to be copied is a vector view, the values
111       of the view will be copied, i.e. the view is not copied.
112
113       \throw A GSL_error is indirectly thrown if memory allocation
114       fails.
115    */
116    vector(const vector& other);
117
118    /**
119       \brief Vector view constructor.
120
121       Create a view of vector \a v, with starting index \a offset,
122       size \a n, and an optional \a stride.
123
124       A vector view can be used as any vector with the difference
125       that changes made to the view will also change the object that
126       is viewed. Also, using the copy constructor will create a new
127       vector object that is a copy of whatever is viewed. If a copy
128       of the view is needed then you should use this constructor to
129       obtain a copy.
130
131       \note If the object viewed by the view goes out of scope or is
132       deleted, the view becomes invalid and the result of further use
133       is undefined.
134
135       \throw GSL_error if a view cannot be set up.
136    */
137    vector(vector& v, size_t offset, size_t n, size_t stride=1);
138
139    /**
140       \brief Vector const view constructor.
141
142       Create a view of vector \a v, with starting index \a offset,
143       size \a n, and an optional \a stride.
144
145       A vector view can be used as any const vector. Using the copy
146       constructor will create a new vector object that is a copy of
147       whatever is viewed. If a copy of the view is needed then you
148       should use this constructor to obtain a copy.
149
150       \note If the object viewed by the view goes out of scope or is
151       deleted, the view becomes invalid and the result of further use
152       is undefined.
153
154       \throw GSL_error if a view cannot be set up.
155    */
156    vector(const vector& v, size_t offset, size_t n, size_t stride=1);
157
158    ///
159    /// Matrix row/column view constructor.
160    ///
161    /// Create a row/column vector view of matrix \a m, pointing at
162    /// row/column \a i. The parameter \a row is used to set whether
163    /// the view should be a row or column view. If \a row is set to
164    /// true, the view will be a row view (default behaviour), and,
165    /// naturally, a column view otherwise.
166    ///
167    /// A vector view can be used as any vector with the difference
168    /// that changes made to the view will also change the object that
169    /// is viewed. Also, using the copy constructor will create a new
170    /// vector object that is a copy of whatever is viewed. If a copy
171    /// of the view is needed then you should use the vector view
172    /// constructor to obtain a copy.
173    ///
174    /// @note If the object viewed by the view goes out of scope or is
175    /// deleted, the view becomes invalid and the result of further
176    /// use is undefined.
177    ///
178    vector(matrix& m, size_t i, bool row=true);
179
180    ///
181    /// Matrix row/column const view constructor.
182    ///
183    /// Create a row/column vector view of matrix \a m, pointing at
184    /// row/column \a i. The parameter \a row is used to set whether
185    /// the view should be a row or column view. If \a row is set to
186    /// true, the view will be a row view (default behaviour), and,
187    /// naturally, a column view otherwise.
188    ///
189    /// A const vector view can be used as any const vector. Using the
190    /// copy constructor will create a new vector object that is a
191    /// copy of whatever is viewed. If a copy of the view is needed
192    /// then you should use the vector view constructor to obtain a
193    /// copy.
194    ///
195    /// @note If the object viewed by the view goes out of scope or is
196    /// deleted, the view becomes invalid and the result of further
197    /// use is undefined.
198    ///
199    vector(const matrix& m, size_t i, bool row=true);
200
201    /**
202       \brief The istream constructor.
203
204       Either elements should be separated with white space characters
205       (default), or elements should be separated by the delimiter \a
206       sep. When delimiter \a sep is used empty elements are stored as
207       NaN's (except that empty lines are ignored). The end of input
208       to the vector is at end of file marker.
209
210       \throw GSL_error if memory allocation fails, IO_error if
211       unexpected input is found in the input stream.
212    */
213    explicit vector(std::istream &, char sep='\0')
214      throw(utility::IO_error, std::exception);
215
216    ///
217    /// The destructor.
218    ///
219    ~vector(void);
220
221    /**
222     */
223    iterator begin(void);
224
225    /**
226     */
227    const_iterator begin(void) const;
228
229    /**
230       \brief Make a copy of \a other.
231
232       This function will make a deep copy of \a other. Memory is
233       resized and view state is changed if needed.
234    */
235    const vector& clone(const vector& other);
236
237    /**
238       \brief This function performs element-wise division, \f$ this_i =
239       this_i/other_i \; \forall i \f$.
240
241       \throw GSL_error if dimensions mis-match.
242    */
243    void div(const vector& other);
244
245    /**
246     */
247    iterator end(void);
248
249    /**
250     */
251    const_iterator end(void) const;
252
253    /**
254       \brief Check whether vectors are equal within a user defined
255       precision, set by \a precision.
256
257       \return True if each element deviates less or equal than \a
258       d. If any vector contain a NaN, false is always returned.
259
260       \see operator== and operator!=
261    */
262    bool equal(const vector&, const double precision=0) const;
263
264    ///
265    /// @return A const pointer to the internal GSL vector,
266    ///
267    const gsl_vector* gsl_vector_p(void) const;
268
269    ///
270    /// @return A pointer to the internal GSL vector,
271    ///
272    gsl_vector* gsl_vector_p(void);
273
274    ///
275    /// Check if the vector object is a view (sub-vector) to another
276    /// vector.
277    ///
278    /// @return True if the object is a view, false othwerwise.
279    ///
280    bool isview(void) const;
281
282    /**
283       \brief This function performs element-wise multiplication, \f$
284       this_i = this_i * other_i \; \forall i \f$.
285
286       \throw GSL_error if dimensions mis-match.
287    */
288    void mul(const vector& other);
289
290    /**
291       @brief Resize vector
292       
293       All elements are set to @a init_value.
294
295       \note Underlying GSL vector is destroyed and a view into this
296       vector becomes invalid.
297    */
298    void resize(size_t, double init_value=0);
299
300    /**
301       \brief Reverse the order of elements in the vector.
302    */
303    void reverse(void);
304
305    ///
306    /// Set all elements to \a value.
307    ///
308    void all(const double& value);
309
310    ///
311    /// @return the number of elements in the vector.
312    ///
313    size_t size(void) const;
314
315    /**
316       \brief Exchange elements \a i and \a j.
317
318       \throw GSL_error if vector lengths differs.
319    */
320    void swap(size_t i, size_t j);
321
322    /**
323       \brief Element access operator.
324
325       \return Reference to element \a i.
326
327       \throw If GSL range checks are enabled in the underlying GSL
328       library a GSL_error exception is thrown if either index is out
329       of range.
330    */
331    double& operator()(size_t i);
332
333    /**
334       \brief Element access operator.
335
336       \return Const reference to element \a i.
337
338       \throw If GSL range checks are enabled in the underlying GSL
339       library a GSL_error exception is thrown if either index is out
340       of range.
341    */
342    const double& operator()(size_t i) const;
343
344    ///
345    /// Element access operator.
346    ///
347    /// @return Reference to element \a i.
348    ///
349    double& operator[](size_t i);
350
351    ///
352    /// Const element access operator.
353    ///
354    /// @return The value of element \a i.
355    ///
356    const double& operator[](size_t i) 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 True if all elements are equal otherwise false.
366
367       \see equal(const vector&, const double precision=0)
368    */
369    bool operator==(const vector&) const;
370
371    /**
372       \brief Comparison operator. Takes linear time.
373
374       Checks are performed with exact matching, i.e., rounding off
375       effects may destroy comparison. Use the equal function for
376       comparing elements within a user defined precision.
377
378       \return False if all elements are equal otherwise true.
379
380       \see equal(const vector&, const double precision=0)
381    */
382    bool operator!=(const vector&) const;
383
384    ///
385    /// @return The dot product.
386    ///
387    double operator*(const vector&) const;
388
389    /**
390       \brief The assignment operator.
391
392       Dimensions of the vectors must match. If the LHS vector is a
393       view, the underlying data will be changed.
394
395       \return A const reference to the resulting vector.
396
397       \see void set(const vector&).
398
399       \throw GSL_error if dimensions mis-match.
400    */
401    const vector& operator=(const vector&);
402
403    /**
404       \brief Addition and assign operator. Vector addition, \f$
405       this_i = this_i + other_i \; \forall i \f$.
406
407       \return A const reference to the resulting vector.
408
409       \throw GSL_error if dimensions mis-match.
410    */
411    const vector& operator+=(const vector&);
412
413    /**
414       \brief Add a constant to a vector, \f$ this_i = this_i + d \;
415       \forall i \f$.
416
417       \return A const reference to the resulting vector.
418    */
419    const vector& operator+=(double d);
420
421    /**
422       \brief Subtract and assign operator. Vector subtraction, \f$
423       this_i = this_i - other_i \; \forall i \f$.
424
425       \return A const reference to the resulting vector.
426
427       \throw GSL_error if dimensions mis-match.
428    */
429    const vector& operator-=(const vector&);
430
431    /**
432       \brief Subtract a constant to a vector, \f$ this_i = this_i - d
433       \; \forall i \f$.
434
435       \return A const reference to the resulting vector.
436    */
437    const vector& operator-=(double d);
438
439    /**
440       \brief Multiply with scalar and assign operator, \f$ this_i =
441       this_i * d \; \forall i \f$.
442
443       \return A const reference to the resulting vector.
444    */
445    const vector& operator*=(double d);
446
447
448  private:
449
450    /**
451       \brief Create a new copy of the internal GSL vector.
452
453       Necessary memory for the new GSL vector is allocated and the
454       caller is responsible for freeing the allocated memory.
455
456       \return A pointer to a copy of the internal GSL vector.
457
458       \throw GSL_error if memory cannot be allocated for the new
459       copy, or if dimensions mis-match.
460    */
461    gsl_vector* create_gsl_vector_copy(void) const;
462
463    /**
464       \brief Clear all dynamically allocated memory.
465
466       Internal utility function.
467    */
468    void delete_allocated_memory(void);
469
470    gsl_vector* v_;
471    gsl_vector_view* view_;
472    const gsl_vector_const_view* view_const_;
473    // proxy_v_ is used to access the proper underlying gsl_vector
474    // in all const member functions. It is not used by design for
475    // non-const vector functions and operators. This is to make sure
476    // that runtime errors occur if a const vector is used in an
477    // inappropriate manner such as on left hand side in assignment
478    // (remember, v_ is null for const vector views).
479    const gsl_vector* proxy_v_;
480  };
481
482  /**
483     \brief Check if all elements of the vector are zero.
484
485     \return True if all elements in the vector is zero, false
486     othwerwise.
487  */
488  bool isnull(const vector&);
489
490  /**
491     \brief Get the maximum value of the vector.
492
493     \return The maximum value of the vector.
494  */
495  double max(const vector&);
496
497  /**
498     \brief Locate the maximum value in the vector.
499
500     \return The index to the maximum value of the vector.
501
502     \note Lower index has precedence.
503  */
504  size_t max_index(const vector&);
505
506  /**
507     \brief Get the minimum value of the vector.
508
509     \return The minimum value of the vector.
510  */
511  double min(const vector&);
512
513  /**
514     \brief Locate the minimum value in the vector.
515
516     \return The index to the minimum value of the vector.
517
518     \note Lower index has precedence.
519  */
520  size_t min_index(const vector&);
521
522  /**
523     \brief Create a vector \a flag indicating NaN's in another vector
524     \a templat.
525
526     The \a flag vector is changed to contain 1's and 0's only. A 1
527     means that the corresponding element in the \a templat vector is
528     valid and a zero means that the corresponding element is a NaN.
529
530     \note Space for vector \a flag is reallocated to fit the size of
531     vector \a templat if sizes mismatch.
532
533     \return True if the \a templat vector contains at least one NaN.
534  */
535  bool nan(const vector& templat, vector& flag);
536
537  /**
538     \brief Transforms a vector to a basis vector.
539
540     All elements are set to zero except the \a i-th element which is
541     set to one.
542  */
543  void set_basis(vector&, size_t i);
544
545  /**
546     Randomly shuffles the elements in vector \a invec
547  */
548  void shuffle(vector& invec); 
549
550  /**
551     Sort the elements in the vector.
552
553     \note Bug in GSL: if the vector contains one ore more NaNs the
554     call never returns.
555  */
556  void sort(vector&);
557
558  /**
559     Create a vector \a sort_index containing the indeces of
560     elements in a another vector \a invec.  The elements of \a
561     sort_index give the index of the vector element which would
562     have been stored in that position if the vector had been sorted
563     in place. The first element of \a sort_index gives the index of the least
564     element in \a invec, and the last element of \a sort_index gives the index of the
565     greatest element in \a invec . The vector \a invec is not changed.
566
567  */
568  void sort_index(std::vector<size_t>& sort_index, const vector& invec);
569
570  /** Similar to sort_index but creates a vector with indices to the \a k
571  smallest elements in \a invec. 
572  */
573  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
574  vector& invec);
575
576  /** Similar to sort_index but creates a vector with indices to the \a k
577  largest elements in \a invec. 
578  */
579  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
580  vector& invec);
581
582 
583
584  /**
585     \brief Calculate the sum of all vector elements.
586
587     \return The sum.
588  */
589  double sum(const vector&);
590
591  /**
592     \brief Swap vector elements by copying.
593
594     The two vectors must have the same length.
595
596     \throw GSL_error if vector lengths differs.
597  */
598  void swap(vector&, vector&);
599
600  /**
601     \brief The output operator for the vector class.
602  */
603  std::ostream& operator<<(std::ostream&, const vector&);
604
605}}} // of namespace utility, yat, and theplu
606
607#endif
Note: See TracBrowser for help on using the repository browser.