From ad17f9570a5f0a35cde44cc206255e889821a5ca Mon Sep 17 00:00:00 2001 From: Arnur Nigmetov Date: Mon, 6 Jun 2016 10:50:37 +0200 Subject: Add actual source from previous repos --- .../include/dnn/geometry/euclidean-fixed.h | 190 +++++++++++++++++++++ 1 file changed, 190 insertions(+) create mode 100644 geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h (limited to 'geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h') diff --git a/geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h b/geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h new file mode 100644 index 0000000..a6ccef7 --- /dev/null +++ b/geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h @@ -0,0 +1,190 @@ +#ifndef DNN_GEOMETRY_EUCLIDEAN_FIXED_H +#define DNN_GEOMETRY_EUCLIDEAN_FIXED_H + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "../parallel/tbb.h" // for dnn::vector<...> + +namespace dnn +{ + // TODO: wrap in another namespace (e.g., euclidean) + + template + struct Point: + boost::addable< Point, + boost::subtractable< Point, + boost::dividable2< Point, Real, + boost::multipliable2< Point, Real > > > >, + public boost::array + { + public: + typedef Real Coordinate; + typedef Real DistanceType; + + + public: + Point(size_t id = 0): id_(id) {} + template + Point(const Point& p, size_t id = 0): + id_(id) { *this = p; } + + static size_t dimension() { return D; } + + // Assign a point of different dimension + template + Point& operator=(const Point& p) { for (size_t i = 0; i < (D < DD ? D : DD); ++i) (*this)[i] = p[i]; if (DD < D) for (size_t i = DD; i < D; ++i) (*this)[i] = 0; return *this; } + + Point& operator+=(const Point& p) { for (size_t i = 0; i < D; ++i) (*this)[i] += p[i]; return *this; } + Point& operator-=(const Point& p) { for (size_t i = 0; i < D; ++i) (*this)[i] -= p[i]; return *this; } + Point& operator/=(Real r) { for (size_t i = 0; i < D; ++i) (*this)[i] /= r; return *this; } + Point& operator*=(Real r) { for (size_t i = 0; i < D; ++i) (*this)[i] *= r; return *this; } + + Real norm2() const { Real n = 0; for (size_t i = 0; i < D; ++i) n += (*this)[i] * (*this)[i]; return n; } + Real max_norm() const + { + Real res = std::fabs((*this)[0]); + for (size_t i = 1; i < D; ++i) + if (std::fabs((*this)[i]) > res) + res = std::fabs((*this)[i]); + return res; + } + + Real l1_norm() const + { + Real res = std::fabs((*this)[0]); + for (size_t i = 1; i < D; ++i) + res += std::fabs((*this)[i]); + return res; + } + + Real lp_norm(const Real p) const + { + assert( !std::isinf(p) ); + if ( p == 1.0 ) + return l1_norm(); + Real res = std::pow(std::fabs((*this)[0]), p); + for (size_t i = 1; i < D; ++i) + res += std::pow(std::fabs((*this)[i]), p); + return std::pow(res, 1.0 / p); + } + + // quick and dirty for now; make generic later + //DistanceType distance(const Point& other) const { return sqrt(sq_distance(other)); } + //DistanceType sq_distance(const Point& other) const { return (other - *this).norm2(); } + + DistanceType distance(const Point& other) const { return (other - *this).max_norm(); } + DistanceType p_distance(const Point& other, const double p) const { return (other - *this).lp_norm(p); } + + size_t id() const { return id_; } + size_t& id() { return id_; } + + private: + friend class boost::serialization::access; + + template + void serialize(Archive& ar, const unsigned int version) { ar & boost::serialization::base_object< boost::array >(*this) & id_; } + + private: + size_t id_; + }; + + template + std::ostream& + operator<<(std::ostream& out, const Point& p) + { out << p[0]; for (size_t i = 1; i < D; ++i) out << " " << p[i]; return out; } + + + template + struct PointTraits; // intentionally undefined; should be specialized for each type + + + template + struct PointTraits< Point > // specialization for dnn::Point + { + typedef Point PointType; + typedef const PointType* PointHandle; + typedef std::vector PointContainer; + + typedef typename PointType::Coordinate Coordinate; + typedef typename PointType::DistanceType DistanceType; + + + static DistanceType + distance(const PointType& p1, const PointType& p2) { if (std::isinf(internal_p)) return p1.distance(p2); else return p1.p_distance(p2, internal_p); } + + static DistanceType + distance(PointHandle p1, PointHandle p2) { return distance(*p1,*p2); } + + static size_t dimension() { return D; } + static Real coordinate(const PointType& p, size_t i) { return p[i]; } + static Real& coordinate(PointType& p, size_t i) { return p[i]; } + static Real coordinate(PointHandle p, size_t i) { return coordinate(*p,i); } + + static size_t id(const PointType& p) { return p.id(); } + static size_t& id(PointType& p) { return p.id(); } + static size_t id(PointHandle p) { return id(*p); } + + static PointHandle + handle(const PointType& p) { return &p; } + static const PointType& + point(PointHandle ph) { return *ph; } + + void swap(PointType& p1, PointType& p2) const { return std::swap(p1, p2); } + + static PointContainer + container(size_t n = 0, const PointType& p = PointType()) { return PointContainer(n, p); } + static typename PointContainer::iterator + iterator(PointContainer& c, PointHandle ph) { return c.begin() + (ph - &c[0]); } + static typename PointContainer::const_iterator + iterator(const PointContainer& c, PointHandle ph) { return c.begin() + (ph - &c[0]); } + + // Internal_p determines which norm will be used in Wasserstein metric (not to + // be confused with wassersteinPower parameter: + // we raise \| p - q \|_{internal_p} to wassersteinPower. + static Real internal_p; + + private: + + friend class boost::serialization::access; + + template + void serialize(Archive& ar, const unsigned int version) {} + + }; + + template + Real PointTraits< Point >::internal_p = std::numeric_limits::infinity(); + + + template + void read_points(const std::string& filename, PointContainer& points) + { + typedef typename boost::range_value::type Point; + typedef typename PointTraits::Coordinate Coordinate; + + std::ifstream in(filename.c_str()); + std::string line; + while(std::getline(in, line)) + { + if (line[0] == '#') continue; // comment line in the file + std::stringstream linestream(line); + Coordinate x; + points.push_back(Point()); + size_t i = 0; + while (linestream >> x) + points.back()[i++] = x; + } + } +} + +#endif -- cgit v1.2.3