source: trunk/yat/utility/VectorBase.cc @ 3780

// $Id: VectorBase.cc 3780 2018-12-06 06:57:02Z peter $

/*
  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
  Copyright (C) 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009, 2011, 2012, 2016, 2017, 2018 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/trac/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "VectorBase.h"
#include "Vector.h"
#include "utility.h"

#include <gsl/gsl_blas.h>
#include <gsl/gsl_sort_vector.h>

#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <numeric>
#include <ostream>
#include <string>
#include <utility>
#include <vector>

namespace theplu {
namespace yat {
namespace utility {


  VectorBase::VectorBase(const gsl_vector* v)
    : const_vec_(v)
  {
  }


  VectorBase::~VectorBase(void)
  {
  }


  VectorBase::const_iterator VectorBase::begin(void) const
  {
    if (const_vec_)
      return const_iterator(&(this->operator()(0)), const_vec_->stride);
    return const_iterator(NULL, 1);
  }


  VectorBase::const_iterator VectorBase::end(void) const
  {
    if (const_vec_)
      return const_iterator(&(this->operator()(0))+const_vec_->stride*size(),
                            const_vec_->stride);
    return const_iterator(NULL, 1);
  }


  bool VectorBase::equal(const VectorBase& other, const double d) const
  {
    if (this==&other)
      return true;
    if (size()!=other.size())
      return false;
    // if gsl error handler disabled, out of bounds index will not
    // abort the program.
    for (size_t i=0; i<size(); ++i)
      if (std::abs( (*this)(i)-other(i) ) > d ||
          std::isnan((*this)(i)) || std::isnan(other(i)) )
        return false;
    return true;
  }


  const gsl_vector* VectorBase::gsl_vector_p(void) const
  {
    return const_vec_;
  }


  size_t VectorBase::size(void) const
  {
    if (!const_vec_)
      return 0;
    return const_vec_->size;
  }


  const double& VectorBase::operator()(size_t i) const
  {
    const double* d=gsl_vector_const_ptr(const_vec_, i);
    if (!d)
      throw utility::GSL_error("VectorBase::operator()",GSL_EINVAL);
    return *d;
  }


  bool VectorBase::operator==(const VectorBase& other) const
  {
    return equal(other);
  }


  bool VectorBase::operator!=(const VectorBase& other) const
  {
    return !equal(other);
  }


  double VectorBase::operator*( const VectorBase &other ) const
  {
    double res = 0.0;;
    gsl_blas_ddot(gsl_vector_p(), other.gsl_vector_p(), &res);
    return res;
  }


  bool isnull(const VectorBase& v)
  {
    return gsl_vector_isnull(v.gsl_vector_p());
  }


  double max(const VectorBase& v)
  {
    return gsl_vector_max(v.gsl_vector_p());
  }


  size_t max_index(const VectorBase& v)
  {
    return gsl_vector_max_index(v.gsl_vector_p());
  }


  double min(const VectorBase& v)
  {
    return gsl_vector_min(v.gsl_vector_p());
  }


  size_t min_index(const VectorBase& v)
  {
    return gsl_vector_min_index(v.gsl_vector_p());
  }


  bool nan(const VectorBase& templat, Vector& flag)
  {
    flag.resize(templat.size());
    return binary_weight(templat.begin(), templat.end(), flag.begin());
  }


  void sort_index(std::vector<size_t>& sort_index, const VectorBase& invec)
  {
    assert(invec.gsl_vector_p());
    gsl_permutation* p = gsl_permutation_alloc(invec.size());
    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
    if (status) {
      gsl_permutation_free(p);
      throw utility::GSL_error("sort_index(vector&,const VectorBase&)",status);
    }
    std::vector<size_t> tmp(p->data,p->data+p->size);
    sort_index.swap(tmp);
    gsl_permutation_free(p);
  }


  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
                            const VectorBase& invec)
  {
    assert(invec.gsl_vector_p());
    assert(k<=invec.size());
    sort_index.resize(k);
    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
  }

  void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
                            const VectorBase& invec)
  {
    assert(invec.gsl_vector_p());
    assert(k<=invec.size());
    sort_index.resize(k);
    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
  }


  double sum(const VectorBase& v)
  {
    return std::accumulate(v.begin(), v.end(), 0.0);
  }


  std::ostream& operator<<(std::ostream& s, const VectorBase& a)
  {
    s.setf(std::ios::dec);
    std::streamsize precision = s.precision();
    s.precision(std::numeric_limits<double>().digits10);
    for (size_t j = 0; j < a.size(); ++j) {
      if (j)
        s << s.fill();
      s << a(j);
    }
    s.precision(precision);
    return s;
  }

  namespace detail
  {
    gsl_vector* create_gsl_vector_copy(const gsl_vector* other)
    {
      if (!other)
        return NULL;
      gsl_vector* vec = create_gsl_vector(other->size);
      if (gsl_vector_memcpy(vec, other))
        throw utility::GSL_error("create_gsl_vector_copy memcpy failed");
      return vec;
    }


    gsl_vector* create_gsl_vector(size_t n)
    {
      if (!n)
        return NULL;
      gsl_vector* vec = gsl_vector_alloc(n);
      if (!vec)
        throw utility::GSL_error("gsl_vector_alloc failed");
      return vec;
    }


    gsl_vector* create_gsl_vector(size_t n, double value)
    {
      if (!n)
        return NULL;
      if (value) {
        gsl_vector* vec = create_gsl_vector(n);
        gsl_vector_set_all(vec, value);
        return vec;
      }
      gsl_vector* vec = gsl_vector_calloc(n);
      if (!vec)
        throw utility::GSL_error("gsl_vector_calloc failed");
      return vec;
    }


    bool overlap(const gsl_vector* a, const gsl_vector* b)
    {
      assert(a);
      assert(b);
      assert(a->size);
      assert(b->size);
      const double* a_first = gsl_vector_const_ptr(a, 0);
      const double* a_last = gsl_vector_const_ptr(a, a->size-1);
      const double* b_first = gsl_vector_const_ptr(b, 0);
      const double* b_last = gsl_vector_const_ptr(b, b->size-1);

      return (a_last >= b_first) && (b_last >= a_first);
    }


    bool serial_overlap(const gsl_vector* a, const gsl_vector* b)
    {
      assert(a);
      assert(b);
      assert(a->size);
      assert(b->size);
      assert(a->size == b->size);
      /*
        ABCDEFGHIJKLMNOPQRSTUVWXYZ
        0ab
        1  b a
        2   b    a
        3    b       a
        4     b          a
        5      b             a
        6       b                a

        a_first in mem pos A
        a_last in mem pos Y
        b_first in mem pos B
        b_last in mem pos H

        Both a[1] and b[3] point to mem pos E

        From this follows that a mem pos is problematic if it
        fullfils the following:
        1) It lies inside [a_first, a_last]
        2) It lies inside [b_first, b_last]

        3) The b-line is below the a-line for that position, which is
        a linear condition, hence no intermediate extreme point, and
        we only need to look at the boundaries.
      */

      const double* a_first = gsl_vector_const_ptr(a, 0);
      const double* a_last = gsl_vector_const_ptr(a, a->size-1);
      const double* b_first = gsl_vector_const_ptr(b, 0);
      const double* b_last = gsl_vector_const_ptr(b, b->size-1);
      if ((a_last >= b_first) && (b_last >= a_first)) {
        const double* first = std::max(a_first, b_first);
        const double* last = std::min(a_last, b_last);
        assert(first <= last);

        // check whether b-line is below a-line in first or last,
        // which is equivalent index(b) > index(a)

        // index in x for a is calculated as:
        // (x - a_first) / a->stride
        // so
        // (x - b_first) * a->stride > (x - a_first) * b->stride)
        return ((first - b_first) * a->stride >
                (first - a_first) * b->stride) ||
          ((last - b_first) * a->stride >
           (last - a_first) * b->stride);
      }
      return false;
    }

  }

}}} // of namespace utility, yat, and thep