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

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

changed euclidean to fulfill design criteria, removes sum_squared_deviation from AveragerPairWeighted?, also changed vector::sort to use stl rather than gsl. refs #248 and fixes #253

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.8 KB
Line 
1#ifndef _theplu_yat_utility_vector_
2#define _theplu_yat_utility_vector_
3
4// $Id: vector.h 899 2007-09-26 21:26:25Z 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#include "Iterator.h"
33
34#include <iostream>
35#include <vector>
36#include <utility>
37
38#include <gsl/gsl_vector.h>
39#include <gsl/gsl_sort_vector.h>
40
41namespace theplu {
42namespace yat {
43namespace utility {
44
45  class matrix;
46
47  /**
48     @brief This is the yat interface to GSL vector.
49
50     For time being 'double' is the only type supported.
51
52     \par File streams:
53     Reading and writing vectors to file streams are of course
54     supported. These are implemented without using GSL functionality,
55     and thus binary read and write to streams are not supported.
56
57     \par Vector views:
58     GSL vector views are supported and these are disguised as
59     ordinary utility::vectors. A support function is added,
60     utility::vector::isview(), that can be used to check if a vector
61     object is a view. Note that view vectors do not own the
62     underlying data, and a view is not valid if the vector/matrix
63     owing the data is deallocated.
64
65     \par
66     Currently there is no restriction on how a vector is used when
67     the vector is a const view into another vector or matrix. To
68     avoid unexpected runtime errors, the programmer must declare
69     const view vectors as 'const' in order to get compile time
70     diagnostics about improper use of a const view vector object. If
71     'const' is not used and the const view vector is used erroneously
72     (such as on the left hand side in assignments), the compiler will
73     not catch the error and runtime error will occur. assert(3) is
74     used to catch the runtime error during development. Example on
75     bad and proper use of const view vectors:
76     @code
77  const vector vm(13,1.0);
78  vector v1(vm,2,4);       // bad code! not const!
79  v1(0)=-123;              // accepted by compiler, runtime abort will occur
80                           // or catched by assert depending on compiler flags
81  const vector v2(vm,2,4); // proper code
82  v2(0)=-123;              // not acceptable for the compiler
83     @endcode
84  */
85
86  class vector
87  {
88  public:
89
90    typedef Iterator<double&, vector> iterator;
91    typedef Iterator<const double, const vector> const_iterator;
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     */
222    iterator begin(void);
223
224    /**
225     */
226    const_iterator begin(void) const;
227
228    /**
229       \brief Make a copy of \a other.
230
231       This function will make a deep copy of \a other. Memory is
232       resized and view state is changed if needed.
233    */
234    const vector& clone(const vector& other);
235
236    /**
237       \brief This function performs element-wise division, \f$ this_i =
238       this_i/other_i \; \forall i \f$.
239
240       \throw GSL_error if dimensions mis-match.
241    */
242    void div(const vector& other);
243
244    /**
245     */
246    iterator end(void);
247
248    /**
249     */
250    const_iterator end(void) const;
251
252    /**
253       \brief Check whether vectors are equal within a user defined
254       precision, set by \a precision.
255
256       \return True if each element deviates less or equal than \a
257       d. If any vector contain a NaN, false is always returned.
258
259       \see operator== and operator!=
260    */
261    bool equal(const vector&, const double precision=0) const;
262
263    ///
264    /// @return A const pointer to the internal GSL vector,
265    ///
266    const gsl_vector* gsl_vector_p(void) const;
267
268    ///
269    /// @return A pointer to the internal GSL vector,
270    ///
271    gsl_vector* gsl_vector_p(void);
272
273    ///
274    /// Check if the vector object is a view (sub-vector) to another
275    /// vector.
276    ///
277    /// @return True if the object is a view, false othwerwise.
278    ///
279    bool isview(void) const;
280
281    /**
282       \brief This function performs element-wise multiplication, \f$
283       this_i = this_i * other_i \; \forall i \f$.
284
285       \throw GSL_error if dimensions mis-match.
286    */
287    void mul(const vector& other);
288
289    /**
290       @brief Resize vector
291       
292       All elements are set to @a init_value.
293
294       \note Underlying GSL vector is destroyed and a view into this
295       vector becomes invalid.
296    */
297    void resize(size_t, double init_value=0);
298
299    /**
300       \brief Reverse the order of elements in the vector.
301    */
302    void reverse(void);
303
304    ///
305    /// Set all elements to \a value.
306    ///
307    void all(const double& value);
308
309    ///
310    /// @return the number of elements in the vector.
311    ///
312    size_t size(void) const;
313
314    /**
315       \brief Exchange elements \a i and \a j.
316
317       \throw GSL_error if vector lengths differs.
318    */
319    void swap(size_t i, size_t j);
320
321    /**
322       \brief Element access operator.
323
324       \return Reference to element \a i.
325
326       \throw If GSL range checks are enabled in the underlying GSL
327       library a GSL_error exception is thrown if either index is out
328       of range.
329    */
330    double& operator()(size_t i);
331
332    /**
333       \brief Element access operator.
334
335       \return Const reference to element \a i.
336
337       \throw If GSL range checks are enabled in the underlying GSL
338       library a GSL_error exception is thrown if either index is out
339       of range.
340    */
341    const double& operator()(size_t i) const;
342
343    ///
344    /// Element access operator.
345    ///
346    /// @return Reference to element \a i.
347    ///
348    double& operator[](size_t i);
349
350    ///
351    /// Const element access operator.
352    ///
353    /// @return The value of element \a i.
354    ///
355    const double& operator[](size_t i) const;
356
357    /**
358       \brief Comparison operator. Takes linear time.
359
360       Checks are performed with exact matching, i.e., rounding off
361       effects may destroy comparison. Use the equal function for
362       comparing elements within a user defined precision.
363
364       \return True if all elements are equal otherwise false.
365
366       \see equal(const vector&, const double precision=0)
367    */
368    bool operator==(const vector&) const;
369
370    /**
371       \brief Comparison operator. Takes linear time.
372
373       Checks are performed with exact matching, i.e., rounding off
374       effects may destroy comparison. Use the equal function for
375       comparing elements within a user defined precision.
376
377       \return False if all elements are equal otherwise true.
378
379       \see equal(const vector&, const double precision=0)
380    */
381    bool operator!=(const vector&) const;
382
383    ///
384    /// @return The dot product.
385    ///
386    double operator*(const vector&) const;
387
388    /**
389       \brief The assignment operator.
390
391       Dimensions of the vectors must match. If the LHS vector is a
392       view, the underlying data will be changed.
393
394       \return A const reference to the resulting vector.
395
396       \see void set(const vector&).
397
398       \throw GSL_error if dimensions mis-match.
399    */
400    const vector& operator=(const vector&);
401
402    /**
403       \brief Addition and assign operator. Vector addition, \f$
404       this_i = this_i + other_i \; \forall i \f$.
405
406       \return A const reference to the resulting vector.
407
408       \throw GSL_error if dimensions mis-match.
409    */
410    const vector& operator+=(const vector&);
411
412    /**
413       \brief Add a constant to a vector, \f$ this_i = this_i + d \;
414       \forall i \f$.
415
416       \return A const reference to the resulting vector.
417    */
418    const vector& operator+=(double d);
419
420    /**
421       \brief Subtract and assign operator. Vector subtraction, \f$
422       this_i = this_i - other_i \; \forall i \f$.
423
424       \return A const reference to the resulting vector.
425
426       \throw GSL_error if dimensions mis-match.
427    */
428    const vector& operator-=(const vector&);
429
430    /**
431       \brief Subtract a constant to a vector, \f$ this_i = this_i - d
432       \; \forall i \f$.
433
434       \return A const reference to the resulting vector.
435    */
436    const vector& operator-=(double d);
437
438    /**
439       \brief Multiply with scalar and assign operator, \f$ this_i =
440       this_i * d \; \forall i \f$.
441
442       \return A const reference to the resulting vector.
443    */
444    const vector& operator*=(double d);
445
446
447  private:
448
449    /**
450       \brief Create a new copy of the internal GSL vector.
451
452       Necessary memory for the new GSL vector is allocated and the
453       caller is responsible for freeing the allocated memory.
454
455       \return A pointer to a copy of the internal GSL vector.
456
457       \throw GSL_error if memory cannot be allocated for the new
458       copy, or if dimensions mis-match.
459    */
460    gsl_vector* create_gsl_vector_copy(void) const;
461
462    /**
463       \brief Clear all dynamically allocated memory.
464
465       Internal utility function.
466    */
467    void delete_allocated_memory(void);
468
469    gsl_vector* v_;
470    gsl_vector_view* view_;
471    const gsl_vector_const_view* view_const_;
472    // proxy_v_ is used to access the proper underlying gsl_vector
473    // in all const member functions. It is not used by design for
474    // non-const vector functions and operators. This is to make sure
475    // that runtime errors occur if a const vector is used in an
476    // inappropriate manner such as on left hand side in assignment
477    // (remember, v_ is null for const vector views).
478    const gsl_vector* proxy_v_;
479  };
480
481  /**
482     \brief Check if all elements of the vector are zero.
483
484     \return True if all elements in the vector is zero, false
485     othwerwise.
486  */
487  bool isnull(const vector&);
488
489  /**
490     \brief Get the maximum value of the vector.
491
492     \return The maximum value of the vector.
493  */
494  double max(const vector&);
495
496  /**
497     \brief Locate the maximum value in the vector.
498
499     \return The index to the maximum value of the vector.
500
501     \note Lower index has precedence.
502  */
503  size_t max_index(const vector&);
504
505  /**
506     \brief Get the minimum value of the vector.
507
508     \return The minimum value of the vector.
509  */
510  double min(const vector&);
511
512  /**
513     \brief Locate the minimum value in the vector.
514
515     \return The index to the minimum value of the vector.
516
517     \note Lower index has precedence.
518  */
519  size_t min_index(const vector&);
520
521  /**
522     \brief Create a vector \a flag indicating NaN's in another vector
523     \a templat.
524
525     The \a flag vector is changed to contain 1's and 0's only. A 1
526     means that the corresponding element in the \a templat vector is
527     valid and a zero means that the corresponding element is a NaN.
528
529     \note Space for vector \a flag is reallocated to fit the size of
530     vector \a templat if sizes mismatch.
531
532     \return True if the \a templat vector contains at least one NaN.
533  */
534  bool nan(const vector& templat, vector& flag);
535
536  /**
537     \brief Transforms a vector to a basis vector.
538
539     All elements are set to zero except the \a i-th element which is
540     set to one.
541  */
542  void set_basis(vector&, size_t i);
543
544  /**
545     Randomly shuffles the elements in vector \a invec
546  */
547  void shuffle(vector& invec); 
548
549  /**
550     Sort the elements in the vector.
551  */
552  void sort(vector&);
553
554  /**
555     Create a vector \a sort_index containing the indeces of
556     elements in a another vector \a invec.  The elements of \a
557     sort_index give the index of the vector element which would
558     have been stored in that position if the vector had been sorted
559     in place. The first element of \a sort_index gives the index of the least
560     element in \a invec, and the last element of \a sort_index gives the index of the
561     greatest element in \a invec . The vector \a invec is not changed.
562
563  */
564  void sort_index(std::vector<size_t>& sort_index, const vector& invec);
565
566  /** Similar to sort_index but creates a vector with indices to the \a k
567  smallest elements in \a invec. 
568  */
569  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
570  vector& invec);
571
572  /** Similar to sort_index but creates a vector with indices to the \a k
573  largest elements in \a invec. 
574  */
575  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
576  vector& invec);
577
578 
579
580  /**
581     \brief Calculate the sum of all vector elements.
582
583     \return The sum.
584  */
585  double sum(const vector&);
586
587  /**
588     \brief Swap vector elements by copying.
589
590     The two vectors must have the same length.
591
592     \throw GSL_error if vector lengths differs.
593  */
594  void swap(vector&, vector&);
595
596  /**
597     \brief The output operator for the vector class.
598  */
599  std::ostream& operator<<(std::ostream&, const vector&);
600
601}}} // of namespace utility, yat, and theplu
602
603#endif
Note: See TracBrowser for help on using the repository browser.