yat/utility/VectorBase.cc

Code
Comments
Other
Rev Date Author Line
22 05 Aug 03 peter 1 // $Id$
12 19 Jun 03 daniel 2
570 05 Apr 06 jari 3 /*
570 05 Apr 06 jari 4   Copyright (C) 2003 Daniel Dalevi, Peter Johansson
2119 12 Dec 09 peter 5   Copyright (C) 2004 Jari Häkkinen, Peter Johansson
2119 12 Dec 09 peter 6   Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
2119 12 Dec 09 peter 7   Copyright (C) 2008 Jari Häkkinen, Peter Johansson
3871 01 Mar 20 peter 8   Copyright (C) 2009, 2011, 2012, 2016, 2017, 2018, 2020 Peter Johansson
570 05 Apr 06 jari 9
1437 25 Aug 08 peter 10   This file is part of the yat library, http://dev.thep.lu.se/trac/yat
570 05 Apr 06 jari 11
675 10 Oct 06 jari 12   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 13   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 14   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 15   License, or (at your option) any later version.
570 05 Apr 06 jari 16
675 10 Oct 06 jari 17   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 18   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 19   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
570 05 Apr 06 jari 20   General Public License for more details.
570 05 Apr 06 jari 21
570 05 Apr 06 jari 22   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 23   along with yat. If not, see <http://www.gnu.org/licenses/>.
570 05 Apr 06 jari 24 */
570 05 Apr 06 jari 25
2881 18 Nov 12 peter 26 #include <config.h>
2881 18 Nov 12 peter 27
1015 01 Feb 08 peter 28 #include "VectorBase.h"
1121 22 Feb 08 peter 29 #include "Vector.h"
680 11 Oct 06 jari 30 #include "utility.h"
125 02 Aug 04 peter 31
2660 21 Nov 11 peter 32 #include <gsl/gsl_blas.h>
1680 29 Dec 08 peter 33 #include <gsl/gsl_sort_vector.h>
1680 29 Dec 08 peter 34
899 26 Sep 07 peter 35 #include <algorithm>
783 05 Mar 07 jari 36 #include <cassert>
789 10 Mar 07 jari 37 #include <cmath>
1920 24 Apr 09 peter 38 #include <limits>
3484 23 Mar 16 peter 39 #include <numeric>
1680 29 Dec 08 peter 40 #include <ostream>
1680 29 Dec 08 peter 41 #include <string>
783 05 Mar 07 jari 42 #include <utility>
22 05 Aug 03 peter 43 #include <vector>
42 26 Feb 04 jari 44
42 26 Feb 04 jari 45 namespace theplu {
680 11 Oct 06 jari 46 namespace yat {
616 31 Aug 06 jari 47 namespace utility {
12 19 Jun 03 daniel 48
12 19 Jun 03 daniel 49
1015 01 Feb 08 peter 50   VectorBase::VectorBase(const gsl_vector* v)
1027 02 Feb 08 peter 51     : const_vec_(v)
686 16 Oct 06 jari 52   {
686 16 Oct 06 jari 53   }
12 19 Jun 03 daniel 54
686 16 Oct 06 jari 55
1015 01 Feb 08 peter 56   VectorBase::~VectorBase(void)
686 16 Oct 06 jari 57   {
686 16 Oct 06 jari 58   }
686 16 Oct 06 jari 59
686 16 Oct 06 jari 60
1015 01 Feb 08 peter 61   VectorBase::const_iterator VectorBase::begin(void) const
686 16 Oct 06 jari 62   {
1038 05 Feb 08 peter 63     if (const_vec_)
1038 05 Feb 08 peter 64       return const_iterator(&(this->operator()(0)), const_vec_->stride);
1038 05 Feb 08 peter 65     return const_iterator(NULL, 1);
686 16 Oct 06 jari 66   }
686 16 Oct 06 jari 67
686 16 Oct 06 jari 68
1015 01 Feb 08 peter 69   VectorBase::const_iterator VectorBase::end(void) const
42 26 Feb 04 jari 70   {
1038 05 Feb 08 peter 71     if (const_vec_)
2660 21 Nov 11 peter 72       return const_iterator(&(this->operator()(0))+const_vec_->stride*size(),
1038 05 Feb 08 peter 73                             const_vec_->stride);
1038 05 Feb 08 peter 74     return const_iterator(NULL, 1);
42 26 Feb 04 jari 75   }
12 19 Jun 03 daniel 76
12 19 Jun 03 daniel 77
1015 01 Feb 08 peter 78   bool VectorBase::equal(const VectorBase& other, const double d) const
42 26 Feb 04 jari 79   {
789 10 Mar 07 jari 80     if (this==&other)
789 10 Mar 07 jari 81       return true;
789 10 Mar 07 jari 82     if (size()!=other.size())
789 10 Mar 07 jari 83       return false;
789 10 Mar 07 jari 84     // if gsl error handler disabled, out of bounds index will not
789 10 Mar 07 jari 85     // abort the program.
789 10 Mar 07 jari 86     for (size_t i=0; i<size(); ++i)
1682 29 Dec 08 peter 87       if (std::abs( (*this)(i)-other(i) ) > d ||
995 30 Nov 07 peter 88           std::isnan((*this)(i)) || std::isnan(other(i)) )
789 10 Mar 07 jari 89         return false;
789 10 Mar 07 jari 90     return true;
789 10 Mar 07 jari 91   }
789 10 Mar 07 jari 92
789 10 Mar 07 jari 93
1015 01 Feb 08 peter 94   const gsl_vector* VectorBase::gsl_vector_p(void) const
714 22 Dec 06 jari 95   {
995 30 Nov 07 peter 96     return const_vec_;
714 22 Dec 06 jari 97   }
714 22 Dec 06 jari 98
714 22 Dec 06 jari 99
1015 01 Feb 08 peter 100   size_t VectorBase::size(void) const
714 22 Dec 06 jari 101   {
1008 01 Feb 08 peter 102     if (!const_vec_)
733 06 Jan 07 peter 103       return 0;
995 30 Nov 07 peter 104     return const_vec_->size;
714 22 Dec 06 jari 105   }
714 22 Dec 06 jari 106
714 22 Dec 06 jari 107
1015 01 Feb 08 peter 108   const double& VectorBase::operator()(size_t i) const
714 22 Dec 06 jari 109   {
995 30 Nov 07 peter 110     const double* d=gsl_vector_const_ptr(const_vec_, i);
755 18 Feb 07 jari 111     if (!d)
1015 01 Feb 08 peter 112       throw utility::GSL_error("VectorBase::operator()",GSL_EINVAL);
755 18 Feb 07 jari 113     return *d;
714 22 Dec 06 jari 114   }
714 22 Dec 06 jari 115
714 22 Dec 06 jari 116
1015 01 Feb 08 peter 117   bool VectorBase::operator==(const VectorBase& other) const
714 22 Dec 06 jari 118   {
789 10 Mar 07 jari 119     return equal(other);
341 07 Jun 05 jari 120   }
341 07 Jun 05 jari 121
341 07 Jun 05 jari 122
1015 01 Feb 08 peter 123   bool VectorBase::operator!=(const VectorBase& other) const
789 10 Mar 07 jari 124   {
789 10 Mar 07 jari 125     return !equal(other);
789 10 Mar 07 jari 126   }
789 10 Mar 07 jari 127
789 10 Mar 07 jari 128
1015 01 Feb 08 peter 129   double VectorBase::operator*( const VectorBase &other ) const
475 22 Dec 05 peter 130   {
475 22 Dec 05 peter 131     double res = 0.0;;
2660 21 Nov 11 peter 132     gsl_blas_ddot(gsl_vector_p(), other.gsl_vector_p(), &res);
475 22 Dec 05 peter 133     return res;
475 22 Dec 05 peter 134   }
475 22 Dec 05 peter 135
475 22 Dec 05 peter 136
1015 01 Feb 08 peter 137   bool isnull(const VectorBase& v)
42 26 Feb 04 jari 138   {
782 05 Mar 07 jari 139     return gsl_vector_isnull(v.gsl_vector_p());
782 05 Mar 07 jari 140   }
782 05 Mar 07 jari 141
782 05 Mar 07 jari 142
1015 01 Feb 08 peter 143   double max(const VectorBase& v)
782 05 Mar 07 jari 144   {
782 05 Mar 07 jari 145     return gsl_vector_max(v.gsl_vector_p());
782 05 Mar 07 jari 146   }
782 05 Mar 07 jari 147
782 05 Mar 07 jari 148
1015 01 Feb 08 peter 149   size_t max_index(const VectorBase& v)
782 05 Mar 07 jari 150   {
782 05 Mar 07 jari 151     return gsl_vector_max_index(v.gsl_vector_p());
782 05 Mar 07 jari 152   }
782 05 Mar 07 jari 153
782 05 Mar 07 jari 154
1015 01 Feb 08 peter 155   double min(const VectorBase& v)
782 05 Mar 07 jari 156   {
782 05 Mar 07 jari 157     return gsl_vector_min(v.gsl_vector_p());
782 05 Mar 07 jari 158   }
782 05 Mar 07 jari 159
782 05 Mar 07 jari 160
1015 01 Feb 08 peter 161   size_t min_index(const VectorBase& v)
782 05 Mar 07 jari 162   {
782 05 Mar 07 jari 163     return gsl_vector_min_index(v.gsl_vector_p());
782 05 Mar 07 jari 164   }
782 05 Mar 07 jari 165
782 05 Mar 07 jari 166
3871 01 Mar 20 peter 167   double norm2_squared(const VectorBase& v)
3871 01 Mar 20 peter 168   {
3871 01 Mar 20 peter 169     return v * v;
3871 01 Mar 20 peter 170   }
3871 01 Mar 20 peter 171
3871 01 Mar 20 peter 172
1120 21 Feb 08 peter 173   bool nan(const VectorBase& templat, Vector& flag)
791 11 Mar 07 jari 174   {
1806 18 Feb 09 peter 175     flag.resize(templat.size());
1807 18 Feb 09 peter 176     return binary_weight(templat.begin(), templat.end(), flag.begin());
791 11 Mar 07 jari 177   }
791 11 Mar 07 jari 178
791 11 Mar 07 jari 179
1015 01 Feb 08 peter 180   void sort_index(std::vector<size_t>& sort_index, const VectorBase& invec)
782 05 Mar 07 jari 181   {
877 20 Sep 07 markus 182     assert(invec.gsl_vector_p());
877 20 Sep 07 markus 183     gsl_permutation* p = gsl_permutation_alloc(invec.size());
877 20 Sep 07 markus 184     int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
877 20 Sep 07 markus 185     if (status) {
877 20 Sep 07 markus 186       gsl_permutation_free(p);
3743 12 Jul 18 peter 187       throw utility::GSL_error("sort_index(vector&,const VectorBase&)",status);
877 20 Sep 07 markus 188     }
1519 21 Sep 08 peter 189     std::vector<size_t> tmp(p->data,p->data+p->size);
1519 21 Sep 08 peter 190     sort_index.swap(tmp);
877 20 Sep 07 markus 191     gsl_permutation_free(p);
877 20 Sep 07 markus 192   }
877 20 Sep 07 markus 193
877 20 Sep 07 markus 194
2687 27 Feb 12 peter 195   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
1015 01 Feb 08 peter 196                             const VectorBase& invec)
879 21 Sep 07 markus 197   {
879 21 Sep 07 markus 198     assert(invec.gsl_vector_p());
879 21 Sep 07 markus 199     assert(k<=invec.size());
879 21 Sep 07 markus 200     sort_index.resize(k);
879 21 Sep 07 markus 201     gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
879 21 Sep 07 markus 202   }
2687 27 Feb 12 peter 203
2687 27 Feb 12 peter 204   void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
1015 01 Feb 08 peter 205                             const VectorBase& invec)
879 21 Sep 07 markus 206   {
879 21 Sep 07 markus 207     assert(invec.gsl_vector_p());
879 21 Sep 07 markus 208     assert(k<=invec.size());
879 21 Sep 07 markus 209     sort_index.resize(k);
879 21 Sep 07 markus 210     gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
879 21 Sep 07 markus 211   }
877 20 Sep 07 markus 212
879 21 Sep 07 markus 213
1015 01 Feb 08 peter 214   double sum(const VectorBase& v)
782 05 Mar 07 jari 215   {
3484 23 Mar 16 peter 216     return std::accumulate(v.begin(), v.end(), 0.0);
782 05 Mar 07 jari 217   }
782 05 Mar 07 jari 218
782 05 Mar 07 jari 219
1015 01 Feb 08 peter 220   std::ostream& operator<<(std::ostream& s, const VectorBase& a)
782 05 Mar 07 jari 221   {
420 02 Dec 05 jari 222     s.setf(std::ios::dec);
1919 24 Apr 09 peter 223     std::streamsize precision = s.precision();
1920 24 Apr 09 peter 224     s.precision(std::numeric_limits<double>().digits10);
42 26 Feb 04 jari 225     for (size_t j = 0; j < a.size(); ++j) {
3780 06 Dec 18 peter 226       if (j)
3780 06 Dec 18 peter 227         s << s.fill();
995 30 Nov 07 peter 228       s << a(j);
42 26 Feb 04 jari 229     }
1919 24 Apr 09 peter 230     s.precision(precision);
42 26 Feb 04 jari 231     return s;
42 26 Feb 04 jari 232   }
42 26 Feb 04 jari 233
3623 09 Feb 17 peter 234   namespace detail
3623 09 Feb 17 peter 235   {
3623 09 Feb 17 peter 236     gsl_vector* create_gsl_vector_copy(const gsl_vector* other)
3623 09 Feb 17 peter 237     {
3623 09 Feb 17 peter 238       if (!other)
3623 09 Feb 17 peter 239         return NULL;
3623 09 Feb 17 peter 240       gsl_vector* vec = create_gsl_vector(other->size);
3623 09 Feb 17 peter 241       if (gsl_vector_memcpy(vec, other))
3623 09 Feb 17 peter 242         throw utility::GSL_error("create_gsl_vector_copy memcpy failed");
3623 09 Feb 17 peter 243       return vec;
3623 09 Feb 17 peter 244     }
3623 09 Feb 17 peter 245
3623 09 Feb 17 peter 246
3623 09 Feb 17 peter 247     gsl_vector* create_gsl_vector(size_t n)
3623 09 Feb 17 peter 248     {
3623 09 Feb 17 peter 249       if (!n)
3623 09 Feb 17 peter 250         return NULL;
3623 09 Feb 17 peter 251       gsl_vector* vec = gsl_vector_alloc(n);
3623 09 Feb 17 peter 252       if (!vec)
3623 09 Feb 17 peter 253         throw utility::GSL_error("gsl_vector_alloc failed");
3623 09 Feb 17 peter 254       return vec;
3623 09 Feb 17 peter 255     }
3623 09 Feb 17 peter 256
3623 09 Feb 17 peter 257
3623 09 Feb 17 peter 258     gsl_vector* create_gsl_vector(size_t n, double value)
3623 09 Feb 17 peter 259     {
3623 09 Feb 17 peter 260       if (!n)
3623 09 Feb 17 peter 261         return NULL;
3623 09 Feb 17 peter 262       if (value) {
3623 09 Feb 17 peter 263         gsl_vector* vec = create_gsl_vector(n);
3623 09 Feb 17 peter 264         gsl_vector_set_all(vec, value);
3623 09 Feb 17 peter 265         return vec;
3623 09 Feb 17 peter 266       }
3623 09 Feb 17 peter 267       gsl_vector* vec = gsl_vector_calloc(n);
3623 09 Feb 17 peter 268       if (!vec)
3623 09 Feb 17 peter 269         throw utility::GSL_error("gsl_vector_calloc failed");
3623 09 Feb 17 peter 270       return vec;
3623 09 Feb 17 peter 271     }
3623 09 Feb 17 peter 272
3623 09 Feb 17 peter 273
3623 09 Feb 17 peter 274     bool overlap(const gsl_vector* a, const gsl_vector* b)
3623 09 Feb 17 peter 275     {
3623 09 Feb 17 peter 276       assert(a);
3623 09 Feb 17 peter 277       assert(b);
3623 09 Feb 17 peter 278       assert(a->size);
3623 09 Feb 17 peter 279       assert(b->size);
3623 09 Feb 17 peter 280       const double* a_first = gsl_vector_const_ptr(a, 0);
3623 09 Feb 17 peter 281       const double* a_last = gsl_vector_const_ptr(a, a->size-1);
3623 09 Feb 17 peter 282       const double* b_first = gsl_vector_const_ptr(b, 0);
3623 09 Feb 17 peter 283       const double* b_last = gsl_vector_const_ptr(b, b->size-1);
3623 09 Feb 17 peter 284
3623 09 Feb 17 peter 285       return (a_last >= b_first) && (b_last >= a_first);
3623 09 Feb 17 peter 286     }
3623 09 Feb 17 peter 287
3623 09 Feb 17 peter 288
3623 09 Feb 17 peter 289     bool serial_overlap(const gsl_vector* a, const gsl_vector* b)
3623 09 Feb 17 peter 290     {
3623 09 Feb 17 peter 291       assert(a);
3623 09 Feb 17 peter 292       assert(b);
3623 09 Feb 17 peter 293       assert(a->size);
3623 09 Feb 17 peter 294       assert(b->size);
3623 09 Feb 17 peter 295       assert(a->size == b->size);
3623 09 Feb 17 peter 296       /*
3623 09 Feb 17 peter 297         ABCDEFGHIJKLMNOPQRSTUVWXYZ
3623 09 Feb 17 peter 298         0ab
3623 09 Feb 17 peter 299         1  b a
3623 09 Feb 17 peter 300         2   b    a
3623 09 Feb 17 peter 301         3    b       a
3623 09 Feb 17 peter 302         4     b          a
3623 09 Feb 17 peter 303         5      b             a
3623 09 Feb 17 peter 304         6       b                a
3623 09 Feb 17 peter 305
3623 09 Feb 17 peter 306         a_first in mem pos A
3623 09 Feb 17 peter 307         a_last in mem pos Y
3623 09 Feb 17 peter 308         b_first in mem pos B
3623 09 Feb 17 peter 309         b_last in mem pos H
3623 09 Feb 17 peter 310
3623 09 Feb 17 peter 311         Both a[1] and b[3] point to mem pos E
3623 09 Feb 17 peter 312
3623 09 Feb 17 peter 313         From this follows that a mem pos is problematic if it
3623 09 Feb 17 peter 314         fullfils the following:
3623 09 Feb 17 peter 315         1) It lies inside [a_first, a_last]
3623 09 Feb 17 peter 316         2) It lies inside [b_first, b_last]
3623 09 Feb 17 peter 317
3623 09 Feb 17 peter 318         3) The b-line is below the a-line for that position, which is
3623 09 Feb 17 peter 319         a linear condition, hence no intermediate extreme point, and
3623 09 Feb 17 peter 320         we only need to look at the boundaries.
3623 09 Feb 17 peter 321       */
3623 09 Feb 17 peter 322
3623 09 Feb 17 peter 323       const double* a_first = gsl_vector_const_ptr(a, 0);
3623 09 Feb 17 peter 324       const double* a_last = gsl_vector_const_ptr(a, a->size-1);
3623 09 Feb 17 peter 325       const double* b_first = gsl_vector_const_ptr(b, 0);
3623 09 Feb 17 peter 326       const double* b_last = gsl_vector_const_ptr(b, b->size-1);
3623 09 Feb 17 peter 327       if ((a_last >= b_first) && (b_last >= a_first)) {
3623 09 Feb 17 peter 328         const double* first = std::max(a_first, b_first);
3623 09 Feb 17 peter 329         const double* last = std::min(a_last, b_last);
3623 09 Feb 17 peter 330         assert(first <= last);
3623 09 Feb 17 peter 331
3623 09 Feb 17 peter 332         // check whether b-line is below a-line in first or last,
3623 09 Feb 17 peter 333         // which is equivalent index(b) > index(a)
3623 09 Feb 17 peter 334
3623 09 Feb 17 peter 335         // index in x for a is calculated as:
3623 09 Feb 17 peter 336         // (x - a_first) / a->stride
3623 09 Feb 17 peter 337         // so
3623 09 Feb 17 peter 338         // (x - b_first) * a->stride > (x - a_first) * b->stride)
3623 09 Feb 17 peter 339         return ((first - b_first) * a->stride >
3623 09 Feb 17 peter 340                 (first - a_first) * b->stride) ||
3623 09 Feb 17 peter 341           ((last - b_first) * a->stride >
3623 09 Feb 17 peter 342            (last - a_first) * b->stride);
3623 09 Feb 17 peter 343       }
3623 09 Feb 17 peter 344       return false;
3623 09 Feb 17 peter 345     }
3623 09 Feb 17 peter 346
3623 09 Feb 17 peter 347   }
3623 09 Feb 17 peter 348
686 16 Oct 06 jari 349 }}} // of namespace utility, yat, and thep

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