yat/utility/VectorBase.h

Code
Comments
Other
Rev Date Author Line
995 30 Nov 07 peter 1 #ifndef _theplu_yat_utility_vector_base_
995 30 Nov 07 peter 2 #define _theplu_yat_utility_vector_base_
616 31 Aug 06 jari 3
22 05 Aug 03 peter 4 // $Id$
22 05 Aug 03 peter 5
570 05 Apr 06 jari 6 /*
570 05 Apr 06 jari 7   Copyright (C) 2003 Daniel Dalevi, Peter Johansson
2119 12 Dec 09 peter 8   Copyright (C) 2004 Jari Häkkinen, Peter Johansson
2119 12 Dec 09 peter 9   Copyright (C) 2005 Jari Häkkinen, Peter Johansson, Markus Ringnér
4359 23 Aug 23 peter 10   Copyright (C) 2006 Jari Häkkinen
2119 12 Dec 09 peter 11   Copyright (C) 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
2119 12 Dec 09 peter 12   Copyright (C) 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 13   Copyright (C) 2009, 2017, 2020 Peter Johansson
570 05 Apr 06 jari 14
1437 25 Aug 08 peter 15   This file is part of the yat library, http://dev.thep.lu.se/trac/yat
570 05 Apr 06 jari 16
675 10 Oct 06 jari 17   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 18   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 19   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 20   License, or (at your option) any later version.
570 05 Apr 06 jari 21
675 10 Oct 06 jari 22   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 23   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 24   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
570 05 Apr 06 jari 25   General Public License for more details.
570 05 Apr 06 jari 26
570 05 Apr 06 jari 27   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 28   along with yat. If not, see <http://www.gnu.org/licenses/>.
570 05 Apr 06 jari 29 */
570 05 Apr 06 jari 30
3605 27 Jan 17 peter 31 #include "BasicVector.h"
3605 27 Jan 17 peter 32 #include "BLAS_level1.h"// include these as they are expected when using Vector
3605 27 Jan 17 peter 33 #include "BLAS_level2.h"// include these as they are expected when using Vector
1038 05 Feb 08 peter 34 #include "StrideIterator.h"
3605 27 Jan 17 peter 35 #include "VectorExpression.h"
420 02 Dec 05 jari 36
1680 29 Dec 08 peter 37 #include <iosfwd>
63 19 Apr 04 peter 38 #include <vector>
1792 11 Feb 09 peter 39 #include <cstddef> // size_t
12 19 Jun 03 daniel 40
12 19 Jun 03 daniel 41 #include <gsl/gsl_vector.h>
12 19 Jun 03 daniel 42
42 26 Feb 04 jari 43 namespace theplu {
680 11 Oct 06 jari 44 namespace yat {
616 31 Aug 06 jari 45 namespace utility {
2687 27 Feb 12 peter 46
1120 21 Feb 08 peter 47   class Vector;
420 02 Dec 05 jari 48
686 16 Oct 06 jari 49   /**
2687 27 Feb 12 peter 50      @brief This is the yat interface to GSL vector.
686 16 Oct 06 jari 51
1046 06 Feb 08 peter 52      This is an interface class for vectors containing the const
1046 06 Feb 08 peter 53      interface. For mutable functionality see VectorMutable.
686 16 Oct 06 jari 54   */
3605 27 Jan 17 peter 55   class VectorBase : public BasicVector<VectorBase>
42 26 Feb 04 jari 56   {
42 26 Feb 04 jari 57   public:
1551 06 Oct 08 peter 58     /**
1551 06 Oct 08 peter 59        value_type is double
1551 06 Oct 08 peter 60
1551 06 Oct 08 peter 61        \since New in yat 0.5
1551 06 Oct 08 peter 62      */
1551 06 Oct 08 peter 63     typedef double value_type;
1551 06 Oct 08 peter 64
1551 06 Oct 08 peter 65     /**
1551 06 Oct 08 peter 66        const_reference type is const double&
1551 06 Oct 08 peter 67
1551 06 Oct 08 peter 68        \since New in yat 0.5
1551 06 Oct 08 peter 69      */
1551 06 Oct 08 peter 70     typedef const double& const_reference;
1551 06 Oct 08 peter 71
1015 01 Feb 08 peter 72     /// \brief VectorBase::const_iterator
1041 06 Feb 08 peter 73     typedef StrideIterator<const double*> const_iterator;
880 21 Sep 07 peter 74
754 17 Feb 07 jari 75     /**
995 30 Nov 07 peter 76        \brief Constructor.
754 17 Feb 07 jari 77     */
1027 02 Feb 08 peter 78     VectorBase(const gsl_vector* v=NULL);
1008 01 Feb 08 peter 79
42 26 Feb 04 jari 80     ///
42 26 Feb 04 jari 81     /// The destructor.
42 26 Feb 04 jari 82     ///
1015 01 Feb 08 peter 83     virtual ~VectorBase(void);
12 19 Jun 03 daniel 84
754 17 Feb 07 jari 85     /**
1015 01 Feb 08 peter 86        \return read-only iterator to start of VectorBase
880 21 Sep 07 peter 87      */
880 21 Sep 07 peter 88     const_iterator begin(void) const;
880 21 Sep 07 peter 89
880 21 Sep 07 peter 90     /**
1015 01 Feb 08 peter 91        \return read-only iterator to end of VectorBase
880 21 Sep 07 peter 92      */
880 21 Sep 07 peter 93     const_iterator end(void) const;
880 21 Sep 07 peter 94
880 21 Sep 07 peter 95     /**
1015 01 Feb 08 peter 96        \brief Check whether VectorBases are equal within a user defined
789 10 Mar 07 jari 97        precision, set by \a precision.
789 10 Mar 07 jari 98
789 10 Mar 07 jari 99        \return True if each element deviates less or equal than \a
1015 01 Feb 08 peter 100        d. If any VectorBase contain a NaN, false is always returned.
789 10 Mar 07 jari 101
789 10 Mar 07 jari 102        \see operator== and operator!=
789 10 Mar 07 jari 103     */
1015 01 Feb 08 peter 104     bool equal(const VectorBase&, const double precision=0) const;
789 10 Mar 07 jari 105
1046 06 Feb 08 peter 106     /**
1046 06 Feb 08 peter 107        \return A const pointer to the internal GSL vector,
1046 06 Feb 08 peter 108     */
714 22 Dec 06 jari 109     const gsl_vector* gsl_vector_p(void) const;
227 01 Feb 05 jari 110
1008 01 Feb 08 peter 111     /**
1008 01 Feb 08 peter 112        Check if the vector object is a view (sub-vector) to another
1008 01 Feb 08 peter 113        vector.
2687 27 Feb 12 peter 114
1008 01 Feb 08 peter 115        \return True if the object is a view, false othwerwise.
1008 01 Feb 08 peter 116      */
2687 27 Feb 12 peter 117     virtual bool isview(void) const=0;
1008 01 Feb 08 peter 118
1046 06 Feb 08 peter 119     /**
1046 06 Feb 08 peter 120        \return number of elements in the VectorBase.
1046 06 Feb 08 peter 121     */
714 22 Dec 06 jari 122     size_t size(void) const;
12 19 Jun 03 daniel 123
755 18 Feb 07 jari 124     /**
755 18 Feb 07 jari 125        \brief Element access operator.
755 18 Feb 07 jari 126
755 18 Feb 07 jari 127        \return Const reference to element \a i.
755 18 Feb 07 jari 128
755 18 Feb 07 jari 129        \throw If GSL range checks are enabled in the underlying GSL
755 18 Feb 07 jari 130        library a GSL_error exception is thrown if either index is out
755 18 Feb 07 jari 131        of range.
755 18 Feb 07 jari 132     */
714 22 Dec 06 jari 133     const double& operator()(size_t i) const;
420 02 Dec 05 jari 134
789 10 Mar 07 jari 135     /**
789 10 Mar 07 jari 136        \brief Comparison operator. Takes linear time.
789 10 Mar 07 jari 137
789 10 Mar 07 jari 138        Checks are performed with exact matching, i.e., rounding off
789 10 Mar 07 jari 139        effects may destroy comparison. Use the equal function for
789 10 Mar 07 jari 140        comparing elements within a user defined precision.
789 10 Mar 07 jari 141
789 10 Mar 07 jari 142        \return True if all elements are equal otherwise false.
789 10 Mar 07 jari 143
1015 01 Feb 08 peter 144        \see equal(const VectorBase&, const double precision=0)
789 10 Mar 07 jari 145     */
1015 01 Feb 08 peter 146     bool operator==(const VectorBase&) const;
276 14 Apr 05 peter 147
789 10 Mar 07 jari 148     /**
789 10 Mar 07 jari 149        \brief Comparison operator. Takes linear time.
789 10 Mar 07 jari 150
789 10 Mar 07 jari 151        Checks are performed with exact matching, i.e., rounding off
789 10 Mar 07 jari 152        effects may destroy comparison. Use the equal function for
789 10 Mar 07 jari 153        comparing elements within a user defined precision.
789 10 Mar 07 jari 154
789 10 Mar 07 jari 155        \return False if all elements are equal otherwise true.
789 10 Mar 07 jari 156
1015 01 Feb 08 peter 157        \see equal(const VectorBase&, const double precision=0)
789 10 Mar 07 jari 158     */
1015 01 Feb 08 peter 159     bool operator!=(const VectorBase&) const;
789 10 Mar 07 jari 160
341 07 Jun 05 jari 161     ///
714 22 Dec 06 jari 162     /// @return The dot product.
714 22 Dec 06 jari 163     ///
1015 01 Feb 08 peter 164     double operator*(const VectorBase&) const;
714 22 Dec 06 jari 165
995 30 Nov 07 peter 166   protected:
1125 22 Feb 08 peter 167     /// pointer to underlying GSL vector
995 30 Nov 07 peter 168     const gsl_vector* const_vec_;
1008 01 Feb 08 peter 169
1008 01 Feb 08 peter 170   private:
1125 22 Feb 08 peter 171     // copy assignment not allowed
1015 01 Feb 08 peter 172     const VectorBase& operator=(const VectorBase&);
420 02 Dec 05 jari 173   };
12 19 Jun 03 daniel 174
782 05 Mar 07 jari 175   /**
1015 01 Feb 08 peter 176      \brief Check if all elements of the VectorBase are zero.
12 19 Jun 03 daniel 177
1015 01 Feb 08 peter 178      \return True if all elements in the VectorBase is zero, false
782 05 Mar 07 jari 179      othwerwise.
1887 31 Mar 09 peter 180
1887 31 Mar 09 peter 181      \relates VectorBase
782 05 Mar 07 jari 182   */
1015 01 Feb 08 peter 183   bool isnull(const VectorBase&);
12 19 Jun 03 daniel 184
782 05 Mar 07 jari 185   /**
1015 01 Feb 08 peter 186      \brief Get the maximum value of the VectorBase.
782 05 Mar 07 jari 187
1015 01 Feb 08 peter 188      \return The maximum value of the VectorBase.
1887 31 Mar 09 peter 189
1887 31 Mar 09 peter 190      \relates VectorBase
782 05 Mar 07 jari 191   */
1015 01 Feb 08 peter 192   double max(const VectorBase&);
782 05 Mar 07 jari 193
782 05 Mar 07 jari 194   /**
1015 01 Feb 08 peter 195      \brief Locate the maximum value in the VectorBase.
782 05 Mar 07 jari 196
1015 01 Feb 08 peter 197      \return The index to the maximum value of the VectorBase.
782 05 Mar 07 jari 198
782 05 Mar 07 jari 199      \note Lower index has precedence.
1887 31 Mar 09 peter 200
1887 31 Mar 09 peter 201      \relates VectorBase
782 05 Mar 07 jari 202   */
1015 01 Feb 08 peter 203   size_t max_index(const VectorBase&);
782 05 Mar 07 jari 204
782 05 Mar 07 jari 205   /**
1015 01 Feb 08 peter 206      \brief Get the minimum value of the VectorBase.
782 05 Mar 07 jari 207
1015 01 Feb 08 peter 208      \return The minimum value of the VectorBase.
1887 31 Mar 09 peter 209
1887 31 Mar 09 peter 210      \relates VectorBase
782 05 Mar 07 jari 211   */
1015 01 Feb 08 peter 212   double min(const VectorBase&);
782 05 Mar 07 jari 213
782 05 Mar 07 jari 214   /**
1015 01 Feb 08 peter 215      \brief Locate the minimum value in the VectorBase.
782 05 Mar 07 jari 216
1015 01 Feb 08 peter 217      \return The index to the minimum value of the VectorBase.
782 05 Mar 07 jari 218
782 05 Mar 07 jari 219      \note Lower index has precedence.
1887 31 Mar 09 peter 220
1887 31 Mar 09 peter 221      \relates VectorBase
782 05 Mar 07 jari 222   */
1015 01 Feb 08 peter 223   size_t min_index(const VectorBase&);
782 05 Mar 07 jari 224
782 05 Mar 07 jari 225   /**
1015 01 Feb 08 peter 226      \brief Create a VectorBase \a flag indicating NaN's in another VectorBase
791 11 Mar 07 jari 227      \a templat.
791 11 Mar 07 jari 228
1015 01 Feb 08 peter 229      The \a flag VectorBase is changed to contain 1's and 0's only. A 1
1015 01 Feb 08 peter 230      means that the corresponding element in the \a templat VectorBase is
791 11 Mar 07 jari 231      valid and a zero means that the corresponding element is a NaN.
791 11 Mar 07 jari 232
1046 06 Feb 08 peter 233      \note Space for vector \a flag is reallocated to fit the size of
1015 01 Feb 08 peter 234      VectorBase \a templat if sizes mismatch.
791 11 Mar 07 jari 235
1015 01 Feb 08 peter 236      \return True if the \a templat VectorBase contains at least one NaN.
1887 31 Mar 09 peter 237
1887 31 Mar 09 peter 238      \relates VectorBase
791 11 Mar 07 jari 239   */
1120 21 Feb 08 peter 240   bool nan(const VectorBase& templat, Vector& flag);
791 11 Mar 07 jari 241
791 11 Mar 07 jari 242   /**
3871 01 Mar 20 peter 243      \return L2 norm squared, or the dot product of \a v with itself
3871 01 Mar 20 peter 244
3871 01 Mar 20 peter 245      \relates VectorBase
3871 01 Mar 20 peter 246
3871 01 Mar 20 peter 247      \since New in yat 0.18
3871 01 Mar 20 peter 248    */
3871 01 Mar 20 peter 249   double norm2_squared(const VectorBase& v);
3871 01 Mar 20 peter 250
3871 01 Mar 20 peter 251
3871 01 Mar 20 peter 252   /**
1008 01 Feb 08 peter 253      Create a vector \a sort_index containing the indeces of
1015 01 Feb 08 peter 254      elements in a another VectorBase \a invec.  The elements of \a
1015 01 Feb 08 peter 255      sort_index give the index of the VectorBase element which would
1015 01 Feb 08 peter 256      have been stored in that position if the VectorBase had been sorted
877 20 Sep 07 markus 257      in place. The first element of \a sort_index gives the index of the least
2687 27 Feb 12 peter 258      element in \a invec, and the last element of \a sort_index gives the
2687 27 Feb 12 peter 259      index of the greatest element in \a invec . The VectorBase \a invec
995 30 Nov 07 peter 260      is not changed.
1887 31 Mar 09 peter 261
1887 31 Mar 09 peter 262      \relatesalso VectorBase
877 20 Sep 07 markus 263   */
1015 01 Feb 08 peter 264   void sort_index(std::vector<size_t>& sort_index, const VectorBase& invec);
877 20 Sep 07 markus 265
2687 27 Feb 12 peter 266   /**
1046 06 Feb 08 peter 267       Similar to sort_index but creates a VectorBase with indices to
1046 06 Feb 08 peter 268       the \a k smallest elements in \a invec.
1887 31 Mar 09 peter 269
1887 31 Mar 09 peter 270      \relatesalso VectorBase
879 21 Sep 07 markus 271   */
2687 27 Feb 12 peter 272   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
1015 01 Feb 08 peter 273                            const VectorBase& invec);
877 20 Sep 07 markus 274
2687 27 Feb 12 peter 275   /**
1046 06 Feb 08 peter 276       Similar to sort_index but creates a VectorBase with indices to
1046 06 Feb 08 peter 277       the \a k largest elements in \a invec.
1887 31 Mar 09 peter 278
1887 31 Mar 09 peter 279      \relatesalso VectorBase
879 21 Sep 07 markus 280   */
2687 27 Feb 12 peter 281   void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
1015 01 Feb 08 peter 282                           const VectorBase& invec);
879 21 Sep 07 markus 283
877 20 Sep 07 markus 284   /**
1015 01 Feb 08 peter 285      \brief Calculate the sum of all VectorBase elements.
782 05 Mar 07 jari 286
782 05 Mar 07 jari 287      \return The sum.
1887 31 Mar 09 peter 288
1887 31 Mar 09 peter 289      \relates VectorBase
782 05 Mar 07 jari 290   */
1015 01 Feb 08 peter 291   double sum(const VectorBase&);
782 05 Mar 07 jari 292
782 05 Mar 07 jari 293   /**
1015 01 Feb 08 peter 294      \brief The output operator for the VectorBase class.
1437 25 Aug 08 peter 295
1437 25 Aug 08 peter 296      Elements in VectorBase \a v are sent to ostream \a s and
1437 25 Aug 08 peter 297      separated with the fill character of stream \a s, s.fill(). If
1437 25 Aug 08 peter 298      you, for example, want to print the VectorBase \a v with the
1437 25 Aug 08 peter 299      elements separated by a ':', you do so by:
2687 27 Feb 12 peter 300      \verbatim
2687 27 Feb 12 peter 301      s << setfill(':') << v;
1437 25 Aug 08 peter 302      \endverbatim
1887 31 Mar 09 peter 303
1887 31 Mar 09 peter 304      \relates VectorBase
782 05 Mar 07 jari 305   */
1437 25 Aug 08 peter 306   std::ostream& operator<<(std::ostream& s, const VectorBase& v);
782 05 Mar 07 jari 307
3623 09 Feb 17 peter 308   /// \cond IGNORE_DOXYGEN
3623 09 Feb 17 peter 309
3623 09 Feb 17 peter 310   // Some convenience functions used in Vector and friends.
3623 09 Feb 17 peter 311   namespace detail {
3623 09 Feb 17 peter 312
3623 09 Feb 17 peter 313     /**
3623 09 Feb 17 peter 314        \brief Create a new gsl vector
3623 09 Feb 17 peter 315
3623 09 Feb 17 peter 316        Necessary memory for the new GSL vector is allocated and the
3623 09 Feb 17 peter 317        caller is responsible for freeing the allocated memory.
3623 09 Feb 17 peter 318
3623 09 Feb 17 peter 319        \return A pointer to created GSL vector.
3623 09 Feb 17 peter 320
3623 09 Feb 17 peter 321        \throw GSL_error if memory cannot be allocated for the new
3623 09 Feb 17 peter 322        GSL vector.
3623 09 Feb 17 peter 323     */
3623 09 Feb 17 peter 324     gsl_vector* create_gsl_vector(size_t n);
3623 09 Feb 17 peter 325
3623 09 Feb 17 peter 326     /**
3623 09 Feb 17 peter 327        same as above but sets all values to init
3623 09 Feb 17 peter 328      */
3623 09 Feb 17 peter 329     gsl_vector* create_gsl_vector(size_t n, double init);
3623 09 Feb 17 peter 330
3623 09 Feb 17 peter 331     /**
3623 09 Feb 17 peter 332        \brief Create a new copy of the internal GSL vector.
3623 09 Feb 17 peter 333
3623 09 Feb 17 peter 334        Necessary memory for the new GSL vector is allocated and the
3623 09 Feb 17 peter 335        caller is responsible for freeing the allocated memory.
3623 09 Feb 17 peter 336
3623 09 Feb 17 peter 337        \return A pointer to a copy of the internal GSL vector.
3623 09 Feb 17 peter 338
3623 09 Feb 17 peter 339        \throw GSL_error if memory cannot be allocated for the new
3623 09 Feb 17 peter 340        copy, or if dimensions mis-match.
3623 09 Feb 17 peter 341     */
3623 09 Feb 17 peter 342     gsl_vector* create_gsl_vector_copy(const gsl_vector*);
3623 09 Feb 17 peter 343
3623 09 Feb 17 peter 344
3623 09 Feb 17 peter 345     // return true if a and b overlap in memory
3623 09 Feb 17 peter 346     bool overlap(const gsl_vector* a, const gsl_vector* b);
3623 09 Feb 17 peter 347
3623 09 Feb 17 peter 348     /*
3623 09 Feb 17 peter 349       Is a stricter condition that overlap(2). Basically it returns
3623 09 Feb 17 peter 350       false if it is safe to write code like:
3623 09 Feb 17 peter 351       for (size_t i=0; i<a->size; ++i)
3623 09 Feb 17 peter 352         a[i] = foo(b[i])
3623 09 Feb 17 peter 353
3623 09 Feb 17 peter 354       in other words, it returns true if there exists i and j, i<j, such that
3623 09 Feb 17 peter 355       gsl_vector_const_ptr(a, j) == gsl_vector_const_ptr(b, i)
3623 09 Feb 17 peter 356
3623 09 Feb 17 peter 357       For speed reasons function might return true in some cases when
3623 09 Feb 17 peter 358       condition above is not true, when both a and b have stride>1 and
3623 09 Feb 17 peter 359       elements sit between each other.
3623 09 Feb 17 peter 360
3623 09 Feb 17 peter 361       Note, a and b must have same size.
3623 09 Feb 17 peter 362      */
3623 09 Feb 17 peter 363     bool serial_overlap(const gsl_vector* a, const gsl_vector* b);
3623 09 Feb 17 peter 364   }
3623 09 Feb 17 peter 365
3623 09 Feb 17 peter 366   /// \endcond
3623 09 Feb 17 peter 367
686 16 Oct 06 jari 368 }}} // of namespace utility, yat, and theplu
12 19 Jun 03 daniel 369
420 02 Dec 05 jari 370 #endif

Generated on Sat Oct 14 02:25:28 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)