source: trunk/yat/utility/VectorBase.cc

Last change on this file was 3871, checked in by Peter, 19 months ago

closes #932

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