Point.h

Go to the documentation of this file.

// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2024 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_POINT_H
#define YODA_POINT_H
 
#include "YODA/AnalysisObject.h"
#include "YODA/Exceptions.h"
#include "YODA/Transformation.h"
#include "YODA/Utils/MathUtils.h"
#include "YODA/Utils/PointUtils.h"
#include "YODA/Utils/Traits.h"
#include "YODA/Utils/sortedvector.h"
#include "YODA/Utils/ndarray.h"
#include <utility>
#include <algorithm>
#include <iostream>
 
 
namespace YODA {
 
 
  class Point {
  public:
 
    using ValuePair = std::pair<double,double>;
 
 
    virtual ~Point() {};
 
 
 
    virtual size_t dim() const = 0;
 
    virtual double val(size_t i) const = 0;
 
    virtual void setVal(const size_t i, const double val) = 0;
 
 
 
 
    //virtual const std::pair<double,double>& errs(size_t i) const = 0;
    virtual void setErr(const size_t i, const double e) = 0;
 
    virtual void setErrs(const size_t i, const double eminus, const double eplus) = 0;
 
    virtual void setErrs(const size_t i, const std::pair<double,double>& e) = 0;
 
    //virtual double errMinus(size_t i) const = 0;
    virtual void setErrMinus(const size_t i, const double eminus) = 0;
 
    //virtual double errPlus(size_t i) const = 0;
    virtual void setErrPlus(const size_t i, const double eplus) = 0;
 
    //virtual double errAvg(size_t i) const = 0;
 
    // /// Get value minus negative error for direction @a i
    // double min(size_t i) const = 0;
    // /// Get value plus positive error for direction @a i
    // double max(size_t i) const = 0;*/
 
 
 
 
    virtual void set(const size_t i, const double val, const double e) = 0;
 
    virtual void set(const size_t i, const double val, const double eminus, const double eplus) = 0;
 
    virtual void set(const size_t i, const double val, const std::pair<double,double>& e) = 0;
 
 
 
 
    // /// Scaling of direction @a i
    virtual void scale(const size_t i, const double scale) = 0;
 
    //template<typename FN>
    //void scale(size_t i, FN f) = 0;
 
 
  };
 
 
 
 
 
  template<size_t N>
  class PointBase : public Point {
  protected:
 
     // Convenient aliases
     using Pair = std::pair<double,double>;
     using ValList = std::initializer_list<double>;
     using PairList = std::initializer_list<Pair>;
 
     // extract the content type of an array Arr
     template<typename Arr>
     using containedType = std::decay_t<decltype(*std::declval<Arr>().begin())>;
 
     // check if content type of an array Arr is Pair
     template<typename Arr>
     using containsPair = typename std::is_same<containedType<Arr>, Pair>;
 
     // succeeds if T is an iterable container
     template<typename T>
     using isIterable = std::enable_if_t<Iterable<T>::value>;
 
     // succeeds if T is an iterable container and its content type is Pair
     template<typename T, typename U>
     using isIterableWithPair = std::enable_if_t<(Iterable<T>::value && Iterable<U>::value && containsPair<U>::value)>;
 
  public:
 
    // Typedefs
    using NdVal = typename Utils::ndarray<double, N>;
    using NdValPair = typename Utils::ndarray<std::pair<double,double>, N>;
    using DataSize = std::integral_constant<size_t, 3*N>;
 
 
 
    // @brief Default constructor
    PointBase() {
      clear();
    }
 
 
    template <typename ValRange = ValList, typename = isIterable<ValRange>>
    PointBase(ValRange&& val) : _val(std::forward<ValRange>(val)) { }
 
 
    template <typename ValRange = ValList,
              typename PairRange = PairList,
              typename = isIterableWithPair<ValRange,PairRange>>
    PointBase(ValRange&& val, PairRange&& errs)
      : _val(std::forward<ValRange>(val)), _errs(std::forward<PairRange>(errs)) { }
 
    template <typename ValRange = ValList, typename = isIterable<ValRange>>
    PointBase(ValRange&& val, ValRange&& errs)
       : _val(std::forward<ValRange>(val)) {
      if (val.size() != N || errs.size() != N)
        throw RangeError("Expected " + std::to_string(N) + " dimensions.");
      size_t i = 0;
      auto it = std::begin(errs);
      const auto& itEnd = std::end(errs);
      for (; it != itEnd; ++it) {
        _errs[i++] = std::make_pair(*it, *it);
      }
    }
 
 
    template <typename ValRange = ValList, typename = isIterable<ValRange>>
    PointBase(ValRange&& val, ValRange&& errsdn, ValRange&& errsup)
       : _val(std::forward<ValRange>(val)) {
      if (val.size() != N || errsdn.size() != N || errsup.size() != N)
        throw RangeError("Expected " + std::to_string(N) + " dimensions.");
      size_t i = 0;
      auto itdn = std::begin(errsdn);
      auto itup = std::begin(errsup);
      const auto& itEnd = std::end(errsdn);
      for (; itdn != itEnd; ) {
        _errs[i++] = std::make_pair(*itdn++, *itup++);
      }
    }
 
    PointBase(const PointBase& p) : _val(p._val), _errs(p._errs) { }
 
    PointBase(PointBase&& p) : _val(std::move(p._val)), _errs(std::move(p._errs)) { }
 
 
    size_t dim() const { return N; }
 
 
    void _renderYODA(std::ostream& os, const int width = 13) const noexcept {
 
      for (size_t i = 0; i < N; ++i) {
        os << std::setw(width) << std::left << _val[i] << "\t"
           << std::setw(width) << std::left << _errs[i].first  << "\t"
           << std::setw(width) << std::left << _errs[i].second << "\t";
      }
      os << "\n";
 
    }
 
 
 
    void clear() {
      for (size_t i = 0; i < N; ++i) {
        _val[i] = 0;
        _errs[i] = {0.,0.};
      }
    }
 
    PointBase& operator = (const PointBase& p) {
      if (this != &p) {
        _val = p._val;
        _errs = p._errs;
      }
      return *this;
    }
 
    PointBase& operator = (PointBase&& p) {
      if (this != &p) {
        _val = std::move(p._val);
        _errs = std::move(p._errs);
      }
      return *this;
    }
 
 
 
 
  public:
 
 
    NdVal& vals() { return _val; }
 
    const NdVal& vals() const { return _val; }
 
    double val(size_t i) const {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      return _val[i];
    }
 
    void setVal(const NdVal& val) {
      _val = val;
    }
 
    void setVal(const size_t i, const double val) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      _val[i] = val;
    }
 
 
 
 
    NdValPair& errs() {
      return _errs;
    }
 
    const NdValPair& errs() const {
      return _errs;
    }
 
    Pair errs(const size_t i) const {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      return _errs[i];
    }
 
    double errMinus(const size_t i) const {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      return _errs[i].first;
    }
 
    double errPlus(const size_t i) const {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      return _errs[i].second;
    }
 
    // Get the average error along axis @a i
    double errAvg(const size_t i) const {
      return 0.5*(errMinus(i) + errPlus(i));
    }
 
    double min(const size_t i) const {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      return _val[i] - _errs[i].first;
    }
 
    double max(const size_t i) const {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      return _val[i] + _errs[i].second;
    }
 
    void setErr(const size_t i, const double e) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      const double err = fabs(e);
      _errs[i] = { err, err};
    }
 
    void setErrs(const size_t i, const double eminus, const double eplus) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      _errs[i] = { eminus, eplus};
    }
 
    void setErrs(const size_t i, const std::pair<double,double>& e) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      _errs[i] = e;
    }
 
    void setErrMinus(const size_t i, const double eminus) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      _errs[i].first = eminus;
    }
 
    void setErrPlus(const size_t i, const double eplus) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      _errs[i].second = eplus;
    }
 
 
 
    void set(const size_t i, const double val, const double e) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      const double err = fabs(e);
      _val[i] = val;
      _errs[i] = {err,err};
    }
 
    void set(const size_t i, const double val, const double eminus, const double eplus) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      _val[i] = val;
      _errs[i].first = eminus;
      _errs[i].second = eplus;
    }
 
    void set(const size_t i, const double val, const std::pair<double,double>& e) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      _val[i] = val;
      _errs[i] = e;
    }
 
 
 
    void scaleVal(const size_t i, const double scale) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      _val[i] *= scale;
    }
 
    void scaleErr(const size_t i, const double scale) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      _errs[i].first  *= scale;
      _errs[i].second *= scale;
    }
 
    void scale(const size_t i, const double scale) {
      scaleVal(i, scale);
      scaleErr(i, scale);
    }
 
    void scale(const NdVal& scales) {
      for (size_t i = 0; i < N; ++i) {
        scale(i, scales[i]);
      }
    }
 
    void scale(const Trf<N>& trf) {
      trf.transform(_val, _errs);
    }
 
    void scale(const size_t i, const Trf<N>& trf) {
      if (i >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      trf.transform(_val[i], _errs[i]);
    }
    //
    void transform(const size_t i, const Trf<N>& trf) {
      scale(i, trf);
    }
 
 
 
 
    std::vector<double> _serializeContent() const noexcept {
      std::vector<double> rtn;
      rtn.reserve(DataSize::value);
      rtn.insert(rtn.end(), _val.begin(),  _val.end());
      for (auto err : _errs) {
        rtn.push_back(std::move(err.first));
        rtn.push_back(std::move(err.second));
      }
      return rtn;
    }
 
    void _deserializeContent(const std::vector<double>& data) {
 
      if (data.size() != DataSize::value)
        throw UserError("Length of serialized data should be "+std::to_string(DataSize::value)+"!");
 
      for (size_t i = 0; i < N; ++i) {
        _val[i] = data[i];
        _errs[i] = { data[N+i], data[2*N+i] };
      }
 
    }
 
    // @}
 
  protected:
 
 
    NdVal _val;
    NdValPair _errs;
 
 
  };
 
 
 
 
  template <size_t N>
  inline bool operator==(const PointBase<N>& a, const PointBase<N>& b) {
    // Compare valitions
    for (size_t i = 0; i < N; ++i) {
      if ( !fuzzyEquals(a.vals()[i],        b.vals()[i]) )        return false;
      if ( !fuzzyEquals(a.errs()[i].first,  b.errs()[i].first))  return false;
      if ( !fuzzyEquals(a.errs()[i].second, b.errs()[i].second)) return false;
    }
    return true;
  }
 
  template <size_t N>
  inline bool operator!=(const PointBase<N>& a, const PointBase<N>& b) {
    return !(a == b);
  }
 
 
  template <size_t N>
  inline bool operator<(const PointBase<N>& a, const PointBase<N>& b) {
    #define LT_IF_NOT_EQ(a,b) { if (!fuzzyEquals(a, b)) return a < b; }
    for (size_t i = 0; i < N; ++i) {
      LT_IF_NOT_EQ(a.vals()[i],        b.vals()[i]);
      LT_IF_NOT_EQ(a.errs()[i].first,  b.errs()[i].first);
      LT_IF_NOT_EQ(a.errs()[i].second, b.errs()[i].second);
    }
    #undef LT_IF_NOT_EQ
    return false;
  }
 
  template <size_t N>
  inline bool operator<=(const PointBase<N>& a, const PointBase<N>& b) {
    if (a == b) return true;
    return a < b;
  }
 
  template <size_t N>
  inline bool operator>(const PointBase<N>& a, const PointBase<N>& b) {
    return !(a <= b);
  }
 
  template <size_t N>
  inline bool operator>=(const PointBase<N>& a, const PointBase<N>& b) {
    return !(a < b);
  }
 
 
  template <size_t N>
  class PointND : public PointBase<N> {
    using BaseT = PointBase<N>;
    using BaseT::BaseT;
  };
 
 
 
  template <>
  class PointND<1> : public PointBase<1>,
                     public XDirectionMixin<PointND<1>> {
  public:
 
    using BaseT = PointBase<1>;
    using BaseT::BaseT;
 
 
    PointND(double x, double ex=0.0)
      : BaseT({x}, {{ex, ex}}) {  }
 
 
    PointND(double x, double exminus, double explus)
      : BaseT({x}, {{exminus, explus}}) {  }
 
 
    PointND(double x, const std::pair<double,double>& ex)
      : BaseT( {x}, {ex}) {  }
 
 
    PointND(const BaseT& other) : BaseT(other) {}
 
    PointND(BaseT&& other) : BaseT(std::move(other)) {}
 
 
 
 
    double val(size_t i=0) const {
      if (i >= 1) throw RangeError("Invalid axis int, must be in range 0..dim-1");
      return _val[i];
    }
 
 
  };
 
 
 
  template <>
  class PointND<2> : public PointBase<2>,
                     public XDirectionMixin<PointND<2>>,
                     public YDirectionMixin<PointND<2>> {
  public:
 
    using BaseT = PointBase<2>;
    using BaseT::BaseT;
 
 
    PointND(double x, double y, double ex=0.0, double ey=0.0)
      : BaseT( {x,y}, {{ex,ex}, {ey,ey}}) {  }
 
 
    PointND(double x, double y,
                 double exminus, double explus,
                 double eyminus, double eyplus)
      : BaseT({x,y}, {{exminus,explus}, {eyminus,eyplus}}) {  }
 
 
    PointND(double x, double y, const std::pair<double,double>& ex, const std::pair<double,double>& ey)
      : BaseT({x,y}, {ex, ey}) {  }
 
 
    PointND(const BaseT& other) : BaseT(other) { }
 
    PointND(BaseT&& other) : BaseT(std::move(other)) { }
 
 
 
    // @name Manipulations
 
    void scaleXY(double scalex, double scaley) {
      scaleX(scalex);
      scaleY(scaley);
    }
 
 
  };
 
 
 
  template <>
  class PointND<3> : public PointBase<3>,
                     public XDirectionMixin<PointND<3>>,
                     public YDirectionMixin<PointND<3>>,
                     public ZDirectionMixin<PointND<3>> {
  public:
 
    using BaseT = PointBase<3>;
    using BaseT::BaseT;
 
 
    PointND(double x, double y, double z, double ex=0.0, double ey=0.0, double ez=0.0)
      : BaseT({x,y,z}, {{ex,ex}, {ey,ey}, {ez,ez}}) { }
 
 
    PointND(double x, double y, double z,
                 double exminus, double explus,
                 double eyminus, double eyplus,
                 double ezminus, double ezplus)
      : BaseT({x,y,z}, {{exminus,explus}, {eyminus,eyplus}, {ezminus,ezplus}}) { }
 
    PointND(double x, double y, double z,
                 const std::pair<double,double>& ex,
                 const std::pair<double,double>& ey,
                 const std::pair<double,double>& ez)
      : BaseT({x,y,z}, {ex,ey,ez}) {  }
 
 
    PointND(const BaseT& other) : BaseT(other) { }
 
    PointND(BaseT&& other) : BaseT(std::move(other)) { }
 
 
    // @name Manipulations
 
    void scaleXYZ(double scalex, double scaley, double scalez) {
      scaleX(scalex);
      scaleY(scaley);
      scaleZ(scalez);
    }
 
 
  };
 
 
  using Point1D = PointND<1>;
  using Point2D = PointND<2>;
  using Point3D = PointND<3>;
  using Point4D = PointND<4>;
 
}
 
#endif