diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/CMakeLists.txt b/fastlib/branches/fastlib-stl/fastlib/tree/CMakeLists.txt index 9564f1ef54..53d45592b8 100644 --- a/fastlib/branches/fastlib-stl/fastlib/tree/CMakeLists.txt +++ b/fastlib/branches/fastlib-stl/fastlib/tree/CMakeLists.txt @@ -5,8 +5,14 @@ cmake_minimum_required(VERSION 2.8) # Do not include test programs here set(SOURCES bounds.h + lmetric.h + dhrectbound.h + dhrectbound_impl.h + dballbound.h + dballbound_impl.h kdtree.h kdtree_impl.h + kdtree_impl.cc spacetree.h statistic.h ) diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/bounds.h b/fastlib/branches/fastlib-stl/fastlib/tree/bounds.h index 0cfd7115c9..461a40173b 100644 --- a/fastlib/branches/fastlib-stl/fastlib/tree/bounds.h +++ b/fastlib/branches/fastlib-stl/fastlib/tree/bounds.h @@ -43,878 +43,8 @@ #ifndef TREE_BOUNDS_H #define TREE_BOUNDS_H -#include "../la/matrix.h" -#include "../la/la.h" +#include "lmetric.h" +#include "dhrectbound.h" +#include "dballbound.h" -#include "../math/math_lib.h" - -/** - * Hyper-rectangle bound for an L-metric. - * - * Template parameter t_pow is the metric to use; use 2 for Euclidean (L2). - */ -template -class DHrectBound { - public: - static const int PREFERRED_POWER = t_pow; - - private: - DRange *bounds_; - index_t dim_; - - OBJECT_TRAVERSAL(DHrectBound) { - OT_OBJ(dim_); - OT_ALLOC(bounds_, dim_); - }; - - public: - /** - * Initializes to specified dimensionality with each dimension the empty - * set. - */ - void Init(index_t dimension) { - //DEBUG_ASSERT_MSG(dim_ == BIG_BAD_NUMBER, "Already initialized"); - bounds_ = mem::Alloc(dimension); - - dim_ = dimension; - Reset(); - } - - /** - * Makes this (uninitialized) box the average of the two arguments, - * i.e. the max and min of each range is the average of the maxes and mins - * of the arguments. - * - * Added by: Bill March, 5/7 - */ - void AverageBoxesInit(const DHrectBound& box1, const DHrectBound& box2) { - - index_t dim = box1.dim(); - DEBUG_ASSERT(dim == box2.dim()); - - Init(dim); - - for (index_t i = 0; i < dim; i++) { - - DRange range; - range = box1.get(i) + box2.get(i); - range *= 0.5; - bounds_[i] = range; - - } - - } // AverageBoxes() - - /** - * Resets all dimensions to the empty set. - */ - void Reset() { - for (index_t i = 0; i < dim_; i++) { - bounds_[i].InitEmptySet(); - } - } - - /** - * Determines if a point is within this bound. - */ - bool Contains(const Vector& point) const { - for (index_t i = 0; i < point.length(); i++) { - if (!bounds_[i].Contains(point[i])) { - return false; - } - } - - return true; - } - - bool Contains(const Vector& point, const Vector& box) const { - const DRange *a = this->bounds_; - for (index_t i = 0; i < point.length(); i++){ - if (a[i].hi > a[i].lo){ - if (!bounds_[i].Contains(point[i])){ - return false; - } - } else if(point[i] > a[i].hi & point[i] < a[i].lo){ - return false; - } - } - return true; - } - - - /** Gets the dimensionality */ - index_t dim() const { - return dim_; - } - - /** - * Gets the range for a particular dimension. - */ - const DRange& get(index_t i) const { - DEBUG_BOUNDS(i, dim_); - return bounds_[i]; - } - - /** - * Calculates the maximum distance within the rectangle - */ - double CalculateMaxDistanceSq() const { - double max_distance=0; - for (index_t i = 0; i < dim_; i++) { - max_distance+=math::Sqr(bounds_[i].width()); - } - return max_distance; - } - - /** Calculates the midpoint of the range */ - void CalculateMidpoint(Vector *centroid) const { - centroid->Init(dim_); - for (index_t i = 0; i < dim_; i++) { - (*centroid)[i] = bounds_[i].mid(); - } - } - /** Calculates the midpoint of the range */ - void CalculateMidpointOverwrite(Vector *centroid) const { - for (index_t i = 0; i < dim_; i++) { - (*centroid)[i] = bounds_[i].mid(); - } - } - - - /** - * Calculates minimum bound-to-point squared distance. - */ - double MinDistanceSq(const double *mpoint) const { - double sum = 0; - const DRange *mbound = bounds_; - - index_t d = dim_; - - do { - double v = *mpoint; - double v1 = mbound->lo - v; - double v2 = v - mbound->hi; - - v = (v1 + fabs(v1)) + (v2 + fabs(v2)); - - mbound++; - mpoint++; - - sum += math::Pow(v); // v is non-negative - } while (--d); - - return math::Pow<2, t_pow>(sum) / 4; - } - - /** - * Calcualtes minimum bound-to-bound squared distance, with - * an offset between their respective coordinate systems. - */ - double MinDistanceSq(const DHrectBound& other, const Vector& offset) const { - double sum = 0; - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - index_t mdim = dim_; - - DEBUG_SAME_SIZE(dim_, other.dim_); - //Add Debug for offset vector - - for (index_t d = 0; d < mdim; d++) { - double v1 = b[d].lo - offset[d] - a[d].hi; - double v2 = a[d].lo + offset[d] - b[d].lo; - - double v = (v1 + fabs(v1)) + (v2 + fabs(v2)); - - sum += math::Pow(v); - } - - return math::Pow<2, t_pow>(sum) / 4; - } - - /** - * Calculates minimum bound-to-point squared distance. - */ - double MinDistanceSq(const Vector& point) const { - DEBUG_SAME_SIZE(point.length(), dim_); - return MinDistanceSq(point.ptr()); - } - - /** - * Calculates minimum bound-to-bound squared distance. - * - * Example: bound1.MinDistanceSq(other) for minimum squared distance. - */ - double MinDistanceSq(const DHrectBound& other) const { - double sum = 0; - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - index_t mdim = dim_; - - DEBUG_SAME_SIZE(dim_, other.dim_); - - for (index_t d = 0; d < mdim; d++) { - double v1 = b[d].lo - a[d].hi; - double v2 = a[d].lo - b[d].hi; - // We invoke the following: - // x + fabs(x) = max(x * 2, 0) - // (x * 2)^2 / 4 = x^2 - double v = (v1 + fabs(v1)) + (v2 + fabs(v2)); - - sum += math::Pow(v); // v is non-negative - } - - return math::Pow<2, t_pow>(sum) / 4; - } - - - /** - * Calculates maximum bound-to-point squared distance. - */ - double MaxDistanceSq(const Vector& point) const { - double sum = 0; - - DEBUG_SAME_SIZE(point.length(), dim_); - - for (index_t d = 0; d < dim_; d++) { - double v = std::max(point[d] - bounds_[d].lo, bounds_[d].hi - point[d]); - sum += math::Pow(v); // v is non-negative - } - - return math::Pow<2, t_pow>(sum); - } - - /** - * Calculates maximum bound-to-point squared distance. - */ - double MaxDistanceSq(const double *point) const { - double sum = 0; - - for (index_t d = 0; d < dim_; d++) { - double v = std::max(point[d] - bounds_[d].lo, bounds_[d].hi - point[d]); - sum += math::Pow(v); // v is non-negative - } - - return math::Pow<2, t_pow>(sum); - } - - /** - * Computes maximum distance. - */ - double MaxDistanceSq(const DHrectBound& other) const { - double sum = 0; - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - - DEBUG_SAME_SIZE(dim_, other.dim_); - - for (index_t d = 0; d < dim_; d++) { - double v = std::max(b[d].hi - a[d].lo, a[d].hi - b[d].lo); - sum += math::PowAbs(v); // v is non-negative - } - - return math::Pow<2, t_pow>(sum); - } - - /** - * Computes maximum distance with offset - */ - double MaxDistanceSq(const DHrectBound& other, const Vector& offset) const { - double sum = 0; - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - - DEBUG_SAME_SIZE(dim_, other.dim_); - - for (index_t d = 0; d < dim_; d++) { - double v = std::max(b[d].hi + offset[d] - a[d].lo, - a[d].hi - offset[d] - b[d].lo); - sum += math::PowAbs(v); // v is non-negative - } - - return math::Pow<2, t_pow>(sum); - } - - /** - * Computes minimum distance between boxes in periodic coordinate system - */ - - double PeriodicMinDistanceSq(const DHrectBound& other, const Vector& box_size) - const { - double sum = 0; - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - - DEBUG_SAME_SIZE(dim_, other.dim_); - - for (index_t d = 0; d < dim_; d++){ - double v = 0, d1, d2, d3; - d1 = (a[d].hi > a[d].lo | b[d].hi > b[d].lo)* - min(b[d].lo - a[d].hi, a[d].lo-b[d].hi); - d2 = (a[d].hi > a[d].lo & b[d].hi > b[d].lo)* - min(b[d].lo - a[d].hi, a[d].lo-b[d].hi + box_size[d]); - d3 = (a[d].hi > a[d].lo & b[d].hi > b[d].lo)* - min(b[d].lo - a[d].hi + box_size[d], a[d].lo-b[d].hi); - v = (d1 + fabs(d1)) + (d2+ fabs(d2)) + (d3 + fabs(d3)); - sum += math::Pow(v); - } - return math::Pow<2, t_pow>(sum) / 4; - } - - - -double PeriodicMinDistanceSq(const Vector& point, const Vector& box_size) - const { - double sum = 0; - const DRange *b = this->bounds_; - - for (index_t d = 0; d < dim_; d++){ - double a = point[d]; - double v = 0, bh; - bh = b[d].hi-b[d].lo; - bh = bh - floor(bh / box_size[d])*box_size[d]; - a = a - b[d].lo; - a = a - floor(a / box_size[d])*box_size[d]; - if (bh > a){ - v = min(a - bh, box_size[d]-a); - } - sum += math::Pow(v); - } - - return math::Pow<2, t_pow>(sum); - } - - - - - /** - * Computes maximum distance between boxes in periodic coordinate system - */ - - double PeriodicMaxDistanceSq(const DHrectBound& other, const Vector& box_size) - const { - double sum = 0; - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - - DEBUG_SAME_SIZE(dim_, other.dim_); - - for (index_t d = 0; d < dim_; d++){ - double v = box_size[d] / 2.0; - double dh, dl; - dh = a[d].hi - b[d].lo; - dh = dh - floor(dh / box_size[d])*box_size[d]; - dl = b[d].hi - a[d].lo; - dl = dl - floor(dl / box_size[d])*box_size[d]; - v = max(min(dh,v), min(dl, v)); - - sum += math::PowAbs(v); - } - return math::Pow<2, t_pow>(sum); - } - - - double PeriodicMaxDistanceSq(const Vector& point, const Vector& box_size) - const { - double sum = 0; - const DRange *a = this->bounds_; - for (index_t d = 0; d < dim_; d++){ - double b = point[d]; - double v = box_size[d] / 2.0; - double ah, al; - ah = a[d].hi - b; - ah = ah - floor(ah / box_size[d])*box_size[d]; - if (ah < v){ - v = ah; - } else { - al = a[d].lo - b; - al = al - floor(al / box_size[d])*box_size[d]; - if (al > v){ - v = 2*v-al; - } - } - sum += math::PowAbs(v); - } - return math::Pow<2, t_pow>(sum); - } - - - double MaxDelta(const DHrectBound& other, double box_width, int dim) - const{ - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - double result = 0.5*box_width; - double temp = b[dim].hi - a[dim].lo; - temp = temp - floor(temp / box_width)*box_width; - if (temp > box_width / 2){ - temp = b[dim].lo - a[dim].hi; - temp = temp - floor(temp / box_width)*box_width; - if (temp > box_width / 2){ - result = b[dim].hi - a[dim].lo; - result = result - floor(temp / box_width + 1)*box_width; - } - } else { - result = temp; - } - return result; - } - - double MinDelta(const DHrectBound& other, double box_width, int dim) - const{ - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - double result = -0.5*box_width; - double temp = b[dim].hi - a[dim].lo; - temp = temp - floor(temp/ box_width)*box_width; - if (temp > box_width / 2){ - temp = b[dim].hi - a[dim].hi; - temp = temp - floor(temp / box_width)*box_width; - if (temp > box_width / 2){ - result = temp - box_width; - } - } else { - temp = b[dim].hi - a[dim].hi; - result = temp - floor(temp / box_width)*box_width; - } - return result; - } - - - - /** - * Calculates minimum and maximum bound-to-bound squared distance. - */ - DRange RangeDistanceSq(const DHrectBound& other) const { - double sum_lo = 0; - double sum_hi = 0; - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - index_t mdim = dim_; - - DEBUG_SAME_SIZE(dim_, other.dim_); - - for (index_t d = 0; d < mdim; d++) { - double v1 = b[d].lo - a[d].hi; - double v2 = a[d].lo - b[d].hi; - // We invoke the following: - // x + fabs(x) = max(x * 2, 0) - // (x * 2)^2 / 4 = x^2 - double v_lo = (v1 + fabs(v1)) + (v2 + fabs(v2)); - double v_hi = -std::min(v1, v2); - - sum_lo += math::Pow(v_lo); // v_lo is non-negative - sum_hi += math::Pow(v_hi); // v_hi is non-negative - } - - return DRange(math::Pow<2, t_pow>(sum_lo) / 4, - math::Pow<2, t_pow>(sum_hi)); - } - - /** - * Calculates minimum and maximum bound-to-point squared distance. - */ - DRange RangeDistanceSq(const Vector& point) const { - double sum_lo = 0; - double sum_hi = 0; - const double *mpoint = point.ptr(); - const DRange *mbound = bounds_; - - DEBUG_SAME_SIZE(point.length(), dim_); - - index_t d = dim_; - do { - double v = *mpoint; - double v1 = mbound->lo - v; - double v2 = v - mbound->hi; - - sum_lo += math::Pow((v1 + fabs(v1)) + (v2 + fabs(v2))); - sum_hi += math::Pow(-std::min(v1, v2)); - - mpoint++; - mbound++; - } while (--d); - - return DRange(math::Pow<2, t_pow>(sum_lo) / 4, - math::Pow<2, t_pow>(sum_hi)); - } - - /** - * Calculates closest-to-their-midpoint bounding box distance, - * i.e. calculates their midpoint and finds the minimum box-to-point - * distance. - * - * Equivalent to: - * - * other.CalcMidpoint(&other_midpoint) - * return MinDistanceSqToPoint(other_midpoint) - * - */ - double MinToMidSq(const DHrectBound& other) const { - double sum = 0; - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - - DEBUG_SAME_SIZE(dim_, other.dim_); - - for (index_t d = 0; d < dim_; d++) { - double v = b->mid(); - double v1 = a->lo - v; - double v2 = v - a->hi; - - v = (v1 + fabs(v1)) + (v2 + fabs(v2)); - - a++; - b++; - - sum += math::Pow(v); // v is non-negative - } - - return math::Pow<2, t_pow>(sum) / 4; - } - - /** - * Computes minimax distance, where the other node is trying to avoid me. - */ - double MinimaxDistanceSq(const DHrectBound& other) const { - double sum = 0; - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - index_t mdim = dim_; - - DEBUG_SAME_SIZE(dim_, other.dim_); - - for (index_t d = 0; d < mdim; d++) { - double v1 = b[d].hi - a[d].hi; - double v2 = a[d].lo - b[d].lo; - double v = std::max(v1, v2); - v = (v + fabs(v)); /* truncate negatives to zero */ - sum += math::Pow(v); // v is non-negative - } - - return math::Pow<2, t_pow>(sum) / 4; - } - - /** - * Calculates midpoint-to-midpoint bounding box distance. - */ - double MidDistanceSq(const DHrectBound& other) const { - double sum = 0; - const DRange *a = this->bounds_; - const DRange *b = other.bounds_; - - DEBUG_SAME_SIZE(dim_, other.dim_); - - for (index_t d = 0; d < dim_; d++) { - sum += math::PowAbs(a[d].hi + a[d].lo - b[d].hi - b[d].lo); - } - - return math::Pow<2, t_pow>(sum) / 4; - } - - /** - * Expands this region to include a new point. - */ - DHrectBound& operator |= (const Vector& vector) { - DEBUG_SAME_SIZE(vector.length(), dim_); - - for (index_t i = 0; i < dim_; i++) { - bounds_[i] |= vector[i]; - } - - return *this; - } - - /** - * Expands this region to encompass another bound. - */ - DHrectBound& operator |= (const DHrectBound& other) { - DEBUG_SAME_SIZE(other.dim_, dim_); - - for (index_t i = 0; i < dim_; i++) { - bounds_[i] |= other.bounds_[i]; - } - - return *this; - } - - - /** - * Expand this bounding box to encompass another point. Done to - * minimize added volume in periodic coordinates. - */ - DHrectBound& Add(const Vector& other, const Vector& size){ - DEBUG_SAME_SIZE(other.length(), dim_); - // Catch case of uninitialized bounds - if (bounds_[0].hi < 0){ - for (index_t i = 0; i < dim_; i++){ - bounds_[i] |= other[i]; - } - } - - for (index_t i= 0; i < dim_; i++){ - double ah, al; - ah = bounds_[i].hi - other[i]; - al = bounds_[i].lo - other[i]; - ah = ah - floor(ah / size[i])*size[i]; - al = al - floor(al / size[i])*size[i]; - if (ah < al){ - if (size[i] - ah < al){ - bounds_[i].hi = other[i]; - } else { - bounds_[i].lo = other[i]; - } - } - } - return *this; - } - - - /** - * Expand this bounding box in periodic coordinates, minimizing added volume. - */ - DHrectBound& Add(const DHrectBound& other, const Vector& size){ - if (bounds_[0].hi < 0){ - for (index_t i = 0; i < dim_; i++){ - bounds_[i] |= other.bounds_[i]; - } - } - - for (index_t i = 0; i < dim_; i++) { - double ah, al, bh, bl; - ah = bounds_[i].hi; - al = bounds_[i].lo; - bh = other.bounds_[i].hi; - bl = other.bounds_[i].lo; - ah = ah - al; - bh = bh - al; - bl = bl - al; - ah = ah - floor(ah / size[i])*size[i]; - bh = bh - floor(bh / size[i])*size[i]; - bl = bl - floor(bl / size[i])*size[i]; - - if (((bh > ah) & (bh < bl | ah > bl )) || - (bh >= bl & bl > ah & bh < ah -bl + size[i])){ - bounds_[i].hi = other.bounds_[i].hi; - } - - if (bl > ah && ((bl > bh) || (bh >= ah -bl + size[i]))){ - bounds_[i].lo = other.bounds_[i].lo; - } - - if (unlikely(ah > bl & bl > bh)){ - bounds_[i].lo = 0; - bounds_[i].hi = size[i]; - } - - - } - return *this; - } - - - -}; - -/** - * An L_p metric for vector spaces. - * - * A generic Metric class should simply compute the distance between - * two points. An LMetric operates for integer powers on Vector spaces. - */ -template -class LMetric { - public: - /** - * Computes the distance metric between two points. - */ - static double Distance(const Vector& a, const Vector& b) { - return math::Pow<1, t_pow>( - la::RawLMetric(a.length(), a.ptr(), b.ptr())); - } - - /** - * Computes the distance metric between two points, raised to a - * particular power. - * - * This might be faster so that you could get, for instance, squared - * L2 distance. - */ - template - static double PowDistance(const Vector& a, const Vector& b) { - return math::Pow( - la::RawLMetric(a.length(), a.ptr(), b.ptr())); - } -}; - -/** - * Ball bound that works in arbitrary metric spaces. - * - * See LMetric for an example metric template parameter. - * - * To initialize this, set the radius with @c set_radius - * and set the point by initializing @c point() directly. - */ -template, typename TPoint = Vector> -class DBallBound { - public: - typedef TPoint Point; - typedef TMetric Metric; - - private: - double radius_; - TPoint center_; - - OBJECT_TRAVERSAL(DBallBound) { - OT_OBJ(radius_); - OT_OBJ(center_); - } - - public: - double radius() const { - return radius_; - } - - void set_radius(double d) { - radius_ = d; - } - - const TPoint& center() const { - return center_; - } - - TPoint& center() { - return center_; - } - - /** - * Determines if a point is within this bound. - */ - bool Contains(const Point& point) const { - return MidDistance(point) <= radius_; - } - - /** - * Gets the center. - * - * Don't really use this directly. This is only here for consistency - * with DHrectBound, so it can plug in more directly if a "centroid" - * is needed. - */ - void CalculateMidpoint(Point *centroid) const { - ot::InitCopy(centroid, center_); - } - - /** - * Calculates minimum bound-to-point squared distance. - */ - double MinDistance(const Point& point) const { - return math::ClampNonNegative(MidDistance(point) - radius_); - } - - double MinDistanceSq(const Point& point) const { - return math::Pow<2, 1>(MinDistance(point)); - } - - /** - * Calculates minimum bound-to-bound squared distance. - */ - double MinDistance(const DBallBound& other) const { - double delta = MidDistance(other.center_) - radius_ - other.radius_; - return math::ClampNonNegative(delta); - } - - double MinDistanceSq(const DBallBound& other) const { - return math::Pow<2, 1>(MinDistance(other)); - } - - /** - * Computes maximum distance. - */ - double MaxDistance(const Point& point) const { - return MidDistance(point) + radius_; - } - - double MaxDistanceSq(const Point& point) const { - return math::Pow<2, 1>(MaxDistance(point)); - } - - /** - * Computes maximum distance. - */ - double MaxDistance(const DBallBound& other) const { - return MidDistance(other.center_) + radius_ + other.radius_; - } - - double MaxDistanceSq(const DBallBound& other) const { - return math::Pow<2, 1>(MaxDistance(other)); - } - - /** - * Calculates minimum and maximum bound-to-bound squared distance. - * - * Example: bound1.MinDistanceSq(other) for minimum squared distance. - */ - DRange RangeDistance(const DBallBound& other) const { - double delta = MidDistance(other.center_); - double sumradius = radius_ + other.radius_; - return DRange( - math::ClampNonNegative(delta - sumradius), - delta + sumradius); - } - - DRange RangeDistanceSq(const DBallBound& other) const { - double delta = MidDistance(other.center_); - double sumradius = radius_ + other.radius_; - return DRange( - math::Pow<2, 1>(math::ClampNonNegative(delta - sumradius)), - math::Pow<2, 1>(delta + sumradius)); - } - - /** - * Calculates closest-to-their-midpoint bounding box distance, - * i.e. calculates their midpoint and finds the minimum box-to-point - * distance. - * - * Equivalent to: - * - * other.CalcMidpoint(&other_midpoint) - * return MinDistanceSqToPoint(other_midpoint) - * - */ - double MinToMid(const DBallBound& other) const { - double delta = MidDistance(other.center_) - radius_; - return math::ClampNonNegative(delta); - } - - double MinToMidSq(const DBallBound& other) const { - return math::Pow<2, 1>(MinToMid(other)); - } - - /** - * Computes minimax distance, where the other node is trying to avoid me. - */ - double MinimaxDistance(const DBallBound& other) const { - double delta = MidDistance(other.center_) + other.radius_ - radius_; - return math::ClampNonNegative(delta); - } - - double MinimaxDistanceSq(const DBallBound& other) const { - return math::Pow<2, 1>(MinimaxDistance(other)); - } - - /** - * Calculates midpoint-to-midpoint bounding box distance. - */ - double MidDistance(const DBallBound& other) const { - return MidDistance(other.center_); - } - - double MidDistanceSq(const DBallBound& other) const { - return math::Pow<2, 1>(MidDistance(other)); - } - - double MidDistance(const Point& point) const { - return Metric::Distance(center_, point); - } -}; - -#endif +#endif // TREE_BOUNDS_H diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/dballbound.h b/fastlib/branches/fastlib-stl/fastlib/tree/dballbound.h new file mode 100644 index 0000000000..f07ec57567 --- /dev/null +++ b/fastlib/branches/fastlib-stl/fastlib/tree/dballbound.h @@ -0,0 +1,172 @@ +/* MLPACK 0.2 + * + * Copyright (c) 2008, 2009 Alexander Gray, + * Garry Boyer, + * Ryan Riegel, + * Nikolaos Vasiloglou, + * Dongryeol Lee, + * Chip Mappus, + * Nishant Mehta, + * Hua Ouyang, + * Parikshit Ram, + * Long Tran, + * Wee Chin Wong + * + * Copyright (c) 2008, 2009 Georgia Institute of Technology + * + * This program 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 2 of the + * License, or (at your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA + * 02110-1301, USA. + */ +/** + * @file tree/dballbound.h + * + * Bounds that are useful for binary space partitioning trees. + * Interface to a ball bound that works in arbitrary metric spaces. + * + * @experimental + */ + +#ifndef TREE_DBALLBOUND_H +#define TREE_DBALLBOUND_H + +#include "../la/matrix.h" +#include "../la/la.h" + +#include "../math/math_lib.h" +#include "lmetric.h" + +#include + +/** + * Ball bound that works in arbitrary metric spaces. + * + * See LMetric for an example metric template parameter. + * + * To initialize this, set the radius with @c set_radius + * and set the point by initializing @c point() directly. + */ +template, typename TPoint = arma::vec> +class DBallBound { + public: + typedef TPoint Point; + typedef TMetric Metric; + + private: + double radius_; + TPoint center_; + + OBJECT_TRAVERSAL(DBallBound) { + OT_OBJ(radius_); + OT_OBJ(center_); + } + + public: + /*** + * Return the radius of the ball bound. + */ + double radius() const { return radius_; } + + /*** + * Set the radius of the bound. + */ + void set_radius(double d) { radius_ = d; } + + /*** + * Return the center point. + */ + const TPoint& center() const { return center_; } + + /*** + * Return the center point. + */ + TPoint& center() { return center_; } + + /** + * Determines if a point is within this bound. + */ + bool Contains(const Point& point) const; + + /** + * Gets the center. + * + * Don't really use this directly. This is only here for consistency + * with DHrectBound, so it can plug in more directly if a "centroid" + * is needed. + */ + void CalculateMidpoint(Point *centroid) const; + + /** + * Calculates minimum bound-to-point squared distance. + */ + double MinDistance(const Point& point) const; + double MinDistanceSq(const Point& point) const; + + /** + * Calculates minimum bound-to-bound squared distance. + */ + double MinDistance(const DBallBound& other) const; + double MinDistanceSq(const DBallBound& other) const; + + /** + * Computes maximum distance. + */ + double MaxDistance(const Point& point) const; + double MaxDistanceSq(const Point& point) const; + + /** + * Computes maximum distance. + */ + double MaxDistance(const DBallBound& other) const; + double MaxDistanceSq(const DBallBound& other) const; + + /** + * Calculates minimum and maximum bound-to-bound squared distance. + * + * Example: bound1.MinDistanceSq(other) for minimum squared distance. + */ + DRange RangeDistance(const DBallBound& other) const; + DRange RangeDistanceSq(const DBallBound& other) const; + + /** + * Calculates closest-to-their-midpoint bounding box distance, + * i.e. calculates their midpoint and finds the minimum box-to-point + * distance. + * + * Equivalent to: + * + * other.CalcMidpoint(&other_midpoint) + * return MinDistanceSqToPoint(other_midpoint) + * + */ + double MinToMid(const DBallBound& other) const; + double MinToMidSq(const DBallBound& other) const; + + /** + * Computes minimax distance, where the other node is trying to avoid me. + */ + double MinimaxDistance(const DBallBound& other) const; + double MinimaxDistanceSq(const DBallBound& other) const; + + /** + * Calculates midpoint-to-midpoint bounding box distance. + */ + double MidDistance(const DBallBound& other) const; + double MidDistanceSq(const DBallBound& other) const; + double MidDistance(const Point& point) const; +}; + +#include "dballbound_impl.h" + +#endif diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/dballbound_impl.h b/fastlib/branches/fastlib-stl/fastlib/tree/dballbound_impl.h new file mode 100644 index 0000000000..c4987d7f93 --- /dev/null +++ b/fastlib/branches/fastlib-stl/fastlib/tree/dballbound_impl.h @@ -0,0 +1,209 @@ +/* MLPACK 0.2 + * + * Copyright (c) 2008, 2009 Alexander Gray, + * Garry Boyer, + * Ryan Riegel, + * Nikolaos Vasiloglou, + * Dongryeol Lee, + * Chip Mappus, + * Nishant Mehta, + * Hua Ouyang, + * Parikshit Ram, + * Long Tran, + * Wee Chin Wong + * + * Copyright (c) 2008, 2009 Georgia Institute of Technology + * + * This program 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 2 of the + * License, or (at your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA + * 02110-1301, USA. + */ +/** + * @file tree/dballbound_impl.h + * + * Bounds that are useful for binary space partitioning trees. + * Implementation of DBallBound ball bound metric policy class. + * + * @experimental + */ + +#ifndef TREE_DBALLBOUND_IMPL_H +#define TREE_DBALLBOUND_IMPL_H + +#include "lmetric.h" + +#include +#include "../base/arma_compat.h" + +/** + * Determines if a point is within the bound. + */ +template +bool DBallBound::Contains(const Point& point) const { + return MidDistance(point) <= radius_; +} + +/** + * Gets the center. + * + * Don't really use this directly. This is only here for consistency + * with DHrectBound, so it can plug in more directly if a "centroid" + * is needed. + */ +template +void DBallBound::CalculateMidpoint(Point *centroid) const { + ot::InitCopy(centroid, center_); +} + +/** + * Calculates minimum bound-to-point squared distance. + */ +template +double DBallBound::MinDistance(const Point& point) const { + return math::ClampNonNegative(MidDistance(point) - radius_); +} + +template +double DBallBound::MinDistanceSq(const Point& point) const { + return math::Pow<2, 1>(MinDistance(point)); +} + +/** + * Calculates minimum bound-to-bound squared distance. + */ +template +double DBallBound::MinDistance(const DBallBound& other) const { + double delta = MidDistance(other.center_) - radius_ - other.radius_; + return math::ClampNonNegative(delta); +} + +template +double DBallBound::MinDistanceSq(const DBallBound& other) const { + return math::Pow<2, 1>(MinDistance(other)); +} + +/** + * Computes maximum distance. + */ +template +double DBallBound::MaxDistance(const Point& point) const { + return MidDistance(point) + radius_; +} + +template +double DBallBound::MaxDistanceSq(const Point& point) const { + return math::Pow<2, 1>(MaxDistance(point)); +} + +/** + * Computes maximum distance. + */ +template +double DBallBound::MaxDistance(const DBallBound& other) const { + return MidDistance(other.center_) + radius_ + other.radius_; +} + +template +double DBallBound::MaxDistanceSq(const DBallBound& other) const { + return math::Pow<2, 1>(MaxDistance(other)); +} + +/** + * Calculates minimum and maximum bound-to-bound squared distance. + * + * Example: bound1.MinDistanceSq(other) for minimum squared distance. + */ +template +DRange DBallBound::RangeDistance(const DBallBound& other) const { + double delta = MidDistance(other.center_); + double sumradius = radius_ + other.radius_; + return DRange( + math::ClampNonNegative(delta - sumradius), + delta + sumradius); +} + +template +DRange DBallBound::RangeDistanceSq(const DBallBound& other) const { + double delta = MidDistance(other.center_); + double sumradius = radius_ + other.radius_; + return DRange( + math::Pow<2, 1>(math::ClampNonNegative(delta - sumradius)), + math::Pow<2, 1>(delta + sumradius)); +} + +/** + * Calculates closest-to-their-midpoint bounding box distance, + * i.e. calculates their midpoint and finds the minimum box-to-point + * distance. + * + * Equivalent to: + * + * other.CalcMidpoint(&other_midpoint) + * return MinDistanceSqToPoint(other_midpoint) + * + */ +template +double DBallBound::MinToMid(const DBallBound& other) const { + double delta = MidDistance(other.center_) - radius_; + return math::ClampNonNegative(delta); +} + +template +double DBallBound::MinToMidSq(const DBallBound& other) const { + return math::Pow<2, 1>(MinToMid(other)); +} + +/** + * Computes minimax distance, where the other node is trying to avoid me. + */ +template +double DBallBound::MinimaxDistance(const DBallBound& other) const { + double delta = MidDistance(other.center_) + other.radius_ - radius_; + return math::ClampNonNegative(delta); +} + +template +double DBallBound::MinimaxDistanceSq(const DBallBound& other) const { + return math::Pow<2, 1>(MinimaxDistance(other)); +} + +/** + * Calculates midpoint-to-midpoint bounding box distance. + */ +template< > +inline double DBallBound, GenVector >::MidDistance(const GenVector& point) const { + arma::vec tmp1; + arma::vec tmp2; + arma_compat::vectorToVec(center_, tmp1); + arma_compat::vectorToVec(point, tmp2); + return Metric::Distance(tmp1, tmp2); +} + +template +double DBallBound::MidDistance(const DBallBound& other) const { + return MidDistance(other.center_); +} + +template +double DBallBound::MidDistanceSq(const DBallBound& other) const { + return math::Pow<2, 1>(MidDistance(other)); +} + +template +double DBallBound::MidDistance(const Point& point) const { + return Metric::Distance(center_, point); +} + +#endif + diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound.h b/fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound.h new file mode 100644 index 0000000000..ceb96292aa --- /dev/null +++ b/fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound.h @@ -0,0 +1,234 @@ +/* MLPACK 0.2 + * + * Copyright (c) 2008, 2009 Alexander Gray, + * Garry Boyer, + * Ryan Riegel, + * Nikolaos Vasiloglou, + * Dongryeol Lee, + * Chip Mappus, + * Nishant Mehta, + * Hua Ouyang, + * Parikshit Ram, + * Long Tran, + * Wee Chin Wong + * + * Copyright (c) 2008, 2009 Georgia Institute of Technology + * + * This program 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 2 of the + * License, or (at your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA + * 02110-1301, USA. + */ +/** + * @file tree/dhrectbound.h + * + * Bounds that are useful for binary space partitioning trees. + * + * This file describes the interface for the DHrectBound policy, which + * implements a hyperrectangle bound. + * + * @experimental + */ + +#ifndef TREE_DHRECTBOUND_H +#define TREE_DHRECTBOUND_H + +#include + +/** + * Hyper-rectangle bound for an L-metric. + * + * Template parameter t_pow is the metric to use; use 2 for Euclidean (L2). + */ +template +class DHrectBound { + public: + static const int PREFERRED_POWER = t_pow; + + private: + DRange *bounds_; + index_t dim_; + +// OBJECT_TRAVERSAL(DHrectBound) { +// OT_OBJ(dim_); + // OT_ALLOC(bounds_, dim_); +// }; + + public: + /** + * Empty constructor. + */ + DHrectBound(); + + /** + * Initializes to specified dimensionality with each dimension the empty + * set. + */ + DHrectBound(index_t dimension); + + /** + * Destructor: clean up memory. + */ + ~DHrectBound(); + + /** + * Makes this (uninitialized) box the average of the two arguments, + * i.e. the max and min of each range is the average of the maxes and mins + * of the arguments. + * + * Added by: Bill March, 5/7 + */ + void AverageBoxesInit(const DHrectBound& box1, const DHrectBound& box2); + + /** + * Resets all dimensions to the empty set. + */ + void Reset(); + + /** + * Sets the dimensionality of the bound. + */ + void SetSize(index_t dim); + + /** + * Determines if a point is within this bound. + */ + bool Contains(const arma::vec& point) const; + + /** Gets the dimensionality */ + index_t dim() const { return dim_; } + + /** + * Gets the range for a particular dimension. + */ + const DRange operator[](index_t i) const; + + /** + * Calculates the maximum distance within the rectangle + */ + double CalculateMaxDistanceSq() const; + + /** Calculates the midpoint of the range */ + void CalculateMidpoint(arma::vec& centroid) const; + + /** + * Calculates minimum bound-to-bound squared distance, with + * an offset between their respective coordinate systems. + */ + double MinDistanceSq(const DHrectBound& other, const arma::vec& offset) const; + + /** + * Calculates minimum bound-to-point squared distance. + */ + double MinDistanceSq(const arma::vec& point) const; + + /** + * Calculates minimum bound-to-bound squared distance. + * + * Example: bound1.MinDistanceSq(other) for minimum squared distance. + */ + double MinDistanceSq(const DHrectBound& other) const; + + /** + * Calculates maximum bound-to-point squared distance. + */ + double MaxDistanceSq(const arma::vec& point) const; + + /** + * Calculates maximum bound-to-point squared distance. + */ + //double MaxDistanceSq(const double *point) const; + + /** + * Computes maximum distance. + */ + double MaxDistanceSq(const DHrectBound& other) const; + + /** + * Computes maximum distance with offset + */ + double MaxDistanceSq(const DHrectBound& other, const arma::vec& offset) const; + + /** + * Computes minimum distance between boxes in periodic coordinate system + */ + double PeriodicMinDistanceSq(const DHrectBound& other, const arma::vec& box_size) const; + double PeriodicMinDistanceSq(const arma::vec& point, const arma::vec& box_size) const; + + /** + * Computes maximum distance between boxes in periodic coordinate system + */ + double PeriodicMaxDistanceSq(const DHrectBound& other, const arma::vec& box_size) const; + double PeriodicMaxDistanceSq(const arma::vec& point, const arma::vec& box_size) const; + + double MaxDelta(const DHrectBound& other, double box_width, int dim) const; + double MinDelta(const DHrectBound& other, double box_width, int dim) const; + + /** + * Calculates minimum and maximum bound-to-bound squared distance. + */ + DRange RangeDistanceSq(const DHrectBound& other) const; + + /** + * Calculates minimum and maximum bound-to-point squared distance. + */ + DRange RangeDistanceSq(const arma::vec& point) const; + + /** + * Calculates closest-to-their-midpoint bounding box distance, + * i.e. calculates their midpoint and finds the minimum box-to-point + * distance. + * + * Equivalent to: + * + * other.CalcMidpoint(&other_midpoint) + * return MinDistanceSqToPoint(other_midpoint) + * + */ + double MinToMidSq(const DHrectBound& other) const; + + /** + * Computes minimax distance, where the other node is trying to avoid me. + */ + double MinimaxDistanceSq(const DHrectBound& other) const; + + /** + * Calculates midpoint-to-midpoint bounding box distance. + */ + double MidDistanceSq(const DHrectBound& other) const; + + /** + * Expands this region to include a new point. + */ + DHrectBound& operator|=(const arma::vec& vector); + + /** + * Expands this region to encompass another bound. + */ + DHrectBound& operator|=(const DHrectBound& other); + + /** + * Expand this bounding box to encompass another point. Done to + * minimize added volume in periodic coordinates. + */ + DHrectBound& Add(const arma::vec& other, const arma::vec& size); + + /** + * Expand this bounding box in periodic coordinates, minimizing added volume. + */ + DHrectBound& Add(const DHrectBound& other, const arma::vec& size); +}; + +#include "dhrectbound_impl.h" + +#endif diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound_impl.h b/fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound_impl.h new file mode 100644 index 0000000000..fbc75affc9 --- /dev/null +++ b/fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound_impl.h @@ -0,0 +1,687 @@ +/* MLPACK 0.2 + * + * Copyright (c) 2008, 2009 Alexander Gray, + * Garry Boyer, + * Ryan Riegel, + * Nikolaos Vasiloglou, + * Dongryeol Lee, + * Chip Mappus, + * Nishant Mehta, + * Hua Ouyang, + * Parikshit Ram, + * Long Tran, + * Wee Chin Wong + * + * Copyright (c) 2008, 2009 Georgia Institute of Technology + * + * This program 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 2 of the + * License, or (at your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA + * 02110-1301, USA. + */ +/** + * @file tree/dhrectbound_impl.h + * + * Implementation of hyper-rectangle bound policy class. + * Template parameter t_pow is the metric to use; use 2 for Euclidian (L2). + * + * @experimental + */ + +#ifndef TREE_DHRECTBOUND_IMPL_H +#define TREE_DHRECTBOUND_IMPL_H + +#include + +using arma::vec; + +/** + * Empty constructor + */ +template +DHrectBound::DHrectBound() { + bounds_ = NULL; + dim_ = 0; +} + +/** + * Initializes to specified dimensionality with each dimension the empty + * set. + */ +template +DHrectBound::DHrectBound(index_t dimension) { + //DEBUG_ASSERT_MSG(dim_ == BIG_BAD_NUMBER, "Already initialized"); + bounds_ = new DRange[dimension]; + + dim_ = dimension; + Reset(); +} + +/** + * Destructor: clean up memory + */ +template +DHrectBound::~DHrectBound() { + if(bounds_) + delete[] bounds_; +} + +/** + * Makes this (uninitialized) box the average of the two arguments, + * i.e. the max and min of each range is the average of the maxes and mins + * of the arguments. + * + * Added by: Bill March, 5/7 + */ +template +void DHrectBound::AverageBoxesInit(const DHrectBound& box1, const DHrectBound& box2) { + + dim_ = box1.dim(); + DEBUG_ASSERT(dim_ == box2.dim()); + + if(bounds_) + delete[] bounds_; + bounds_ = new DRange[dim_]; + + for (index_t i = 0; i < dim_; i++) { + DRange range; + range = box1.get(i) + box2.get(i); + range *= 0.5; + bounds_[i] = range; + } +} + +/** + * Resets all dimensions to the empty set. + */ +template +void DHrectBound::Reset() { + for (index_t i = 0; i < dim_; i++) { + bounds_[i].InitEmptySet(); + } +} + +/** + * Sets the dimensionality. + */ +template +void DHrectBound::SetSize(index_t dim) { + if(bounds_) + delete[] bounds_; + + bounds_ = new DRange[dim]; + dim_ = dim; + Reset(); +} + +/** + * Determines if a point is within this bound. + */ +template +bool DHrectBound::Contains(const vec& point) const { + for (index_t i = 0; i < point.n_elem; i++) { + if (!bounds_[i].Contains(point(i))) { + return false; + } + } + + return true; +} + +/** + * Gets the range for a particular dimension. + */ +template +const DRange DHrectBound::operator[](index_t i) const { + return bounds_[i]; +} + +/** + * Calculates the maximum distance within the rectangle + */ +template +double DHrectBound::CalculateMaxDistanceSq() const { + double max_distance = 0; + for (index_t i = 0; i < dim_; i++) + max_distance += pow(bounds_[i].width(), 2); + + return max_distance; +} + +/** Calculates the midpoint of the range */ +template +void DHrectBound::CalculateMidpoint(vec& centroid) const { + // set size correctly if necessary + if(!(centroid.n_elem == dim_)) + centroid.set_size(dim_); + + for(index_t i = 0; i < dim_; i++) { + centroid(i) = bounds_[i].mid(); + } +} + +/** + * Calcualtes minimum bound-to-bound squared distance, with + * an offset between their respective coordinate systems. + */ +template +double DHrectBound::MinDistanceSq(const DHrectBound& other, const vec& offset) const { + double sum = 0; + + DEBUG_SAME_SIZE(dim_, other.dim_); + //Add Debug for offset vector + + for(index_t d = 0; d < dim_; d++) { + double v1 = other.bounds_[d].lo - offset[d] - bounds_[d].hi; + double v2 = bounds_[d].lo + offset[d] - other.bounds_[d].lo; + + double v = (v1 + fabs(v1)) + (v2 + fabs(v2)); + + sum += pow(v, (double) t_pow); + } + + return pow(sum, 2.0 / (double) t_pow) / 4.0; +} + +/** + * Calculates minimum bound-to-point squared distance. + */ +template +double DHrectBound::MinDistanceSq(const vec& point) const { + DEBUG_SAME_SIZE(point.n_elem, dim_); + + double sum = 0; + double lower, higher; + for(index_t d = 0; d < dim_; d++) { + lower = bounds_[d].lo - point[d]; // negative if point[d] > bounds_[d] + higher = point[d] - bounds_[d].hi; // negative if point[d] < bounds_[d] + + // since only one of 'lower' or 'higher' is negative, if we add each's + // absolute value to itself and then sum those two, our result is the + // nonnegative half of the equation times two; then we raise to power t_pow + sum += pow((lower + fabs(lower)) + (higher + fabs(higher)), (double) t_pow); + } + // now take the t_pow'th root (but make sure our result is squared); then + // divide by four to cancel out the constant of 2 (which has been squared now) + // that was introduced earlier + return pow(sum, 2.0 / (double) t_pow) / 4.0; +} + +/** + * Calculates minimum bound-to-bound squared distance. + * + * Example: bound1.MinDistanceSq(other) for minimum squared distance. + */ +template +double DHrectBound::MinDistanceSq(const DHrectBound& other) const { + double sum = 0; + index_t mdim = dim_; + + DEBUG_SAME_SIZE(dim_, other.dim_); + + for (index_t d = 0; d < mdim; d++) { + double v1 = other.bounds_[d].lo - bounds_[d].hi; + double v2 = bounds_[d].lo - other.bounds_[d].hi; + // We invoke the following: + // x + fabs(x) = max(x * 2, 0) + // (x * 2)^2 / 4 = x^2 + double v = (v1 + fabs(v1)) + (v2 + fabs(v2)); + + sum += pow(v, (double) t_pow); + } + + return pow(sum, 2.0 / (double) t_pow) / 4.0; +} + +/** + * Calculates maximum bound-to-point squared distance. + */ +template +double DHrectBound::MaxDistanceSq(const vec& point) const { + double sum = 0; + + DEBUG_SAME_SIZE(point.n_elem, dim_); + + double v; + for (index_t d = 0; d < dim_; d++) { + double v = fabs(std::max( + point[d] - bounds_[d].lo, + bounds_[d].hi - point[d])); + sum += pow(v, (double) t_pow); + } + + return pow(sum, 2.0 / (double) t_pow); +} + +/** + * Calculates maximum bound-to-point squared distance. + */ +/* +template +double DHrectBound::MaxDistanceSq(const double *point) const { + double sum = 0; + + for (index_t d = 0; d < dim_; d++) { + double v = std::max(point[d] - bounds_[d].lo, bounds_[d].hi - point[d]); + sum += math::Pow(v); // v is non-negative + } + + return math::Pow<2, t_pow>(sum); +} +*/ + +/** + * Computes maximum distance. + */ +template +double DHrectBound::MaxDistanceSq(const DHrectBound& other) const { + double sum = 0; + + DEBUG_SAME_SIZE(dim_, other.dim_); + + double v; + for(index_t d = 0; d < dim_; d++) { + v = fabs(std::max( + other.bounds_[d].hi - bounds_[d].lo, + bounds_[d].hi - other.bounds_[d].lo)); + sum += pow(v, (double) t_pow); // v is non-negative + } + + return pow(sum, 2.0 / (double) t_pow); +} + +/** + * Computes maximum distance with offset + */ +template +double DHrectBound::MaxDistanceSq(const DHrectBound& other, const vec& offset) const { + double sum = 0; + + DEBUG_SAME_SIZE(dim_, other.dim_); + + for (index_t d = 0; d < dim_; d++) { + double v = fabs(std::max( + other.bounds_[d].hi + offset[d] - bounds_[d].lo, + bounds_[d].hi - offset[d] - other.bounds_[d].lo)); + sum += pow(v, (double) t_pow); // v is non-negative + } + + return pow(sum, 2.0 / (double) t_pow); +} + +/** + * Computes minimum distance between boxes in periodic coordinate system + */ +template +double DHrectBound::PeriodicMinDistanceSq(const DHrectBound& other, const vec& box_size) const { + double sum = 0; + + DEBUG_SAME_SIZE(dim_, other.dim_); + + for (index_t d = 0; d < dim_; d++){ + double v = 0, d1, d2, d3; + d1 = (bounds_[d].hi > bounds_[d].lo | other.bounds_[d].hi > other.bounds_[d].lo) * + min(other.bounds_[d].lo - bounds_[d].hi, bounds_[d].lo - other.bounds_[d].hi); + d2 = (bounds_[d].hi > bounds_[d].lo & other.bounds_[d].hi > other.bounds_[d].lo) * + min(other.bounds_[d].lo - bounds_[d].hi, bounds_[d].lo - other.bounds_[d].hi + box_size[d]); + d3 = (bounds_[d].hi > bounds_[d].lo & other.bounds_[d].hi > other.bounds_[d].lo) * + min(other.bounds_[d].lo - bounds_[d].hi + box_size[d], bounds_[d].lo - other.bounds_[d].hi); + v = (d1 + fabs(d1)) + (d2 + fabs(d2)) + (d3 + fabs(d3)); + sum += pow(v, (double) t_pow); + } + return pow(sum, 2.0 / (double) t_pow) / 4.0; +} + +template +double DHrectBound::PeriodicMinDistanceSq(const vec& point, const vec& box_size) const { + double sum = 0; + + for (index_t d = 0; d < dim_; d++){ + double a = point[d]; + double v = 0, bh; + bh = bounds_[d].hi - bounds_[d].lo; + bh = bh - floor(bh / box_size[d]) * box_size[d]; + a = a - bounds_[d].lo; + a = a - floor(a / box_size[d]) * box_size[d]; + if (bh > a) + v = min(a - bh, box_size[d]-a); + + sum += pow(v, (double) t_pow); + } + + return pow(sum, 2.0 / (double) t_pow); +} + +/** + * Computes maximum distance between boxes in periodic coordinate system + */ +template +double DHrectBound::PeriodicMaxDistanceSq(const DHrectBound& other, const vec& box_size) const { + double sum = 0; + + DEBUG_SAME_SIZE(dim_, other.dim_); + + for (index_t d = 0; d < dim_; d++){ + double v = box_size[d] / 2.0; + double dh, dl; + dh = bounds_[d].hi - other.bounds_[d].lo; + dh = dh - floor(dh / box_size[d]) * box_size[d]; + dl = other.bounds_[d].hi - bounds_[d].lo; + dl = dl - floor(dl / box_size[d]) * box_size[d]; + v = fabs(max(min(dh, v), min(dl, v))); + + sum += pow(v, (double) t_pow); + } + return pow(sum, 2.0 / (double) t_pow); +} + +template +double DHrectBound::PeriodicMaxDistanceSq(const vec& point, const vec& box_size) const { + double sum = 0; + + for (index_t d = 0; d < dim_; d++) { + double b = point[d]; + double v = box_size[d] / 2.0; + double ah, al; + ah = bounds_[d].hi - b; + ah = ah - floor(ah / box_size[d]) * box_size[d]; + if (ah < v) { + v = ah; + } else { + al = bounds_[d].lo - b; + al = al - floor(al / box_size[d]) * box_size[d]; + if (al > v) { + v = (2 * v) - al; + } + } + sum += pow(fabs(v), (double) t_pow); + } + return pow(sum, 2.0 / (double) t_pow); +} + + +template +double DHrectBound::MaxDelta(const DHrectBound& other, double box_width, int dim) const { + double result = 0.5 * box_width; + double temp = other.bounds_[dim].hi - bounds_[dim].lo; + temp = temp - floor(temp / box_width) * box_width; + if (temp > box_width / 2) { + temp = other.bounds_[dim].lo - bounds_[dim].hi; + temp = temp - floor(temp / box_width) * box_width; + if (temp > box_width / 2) { + result = other.bounds_[dim].hi - bounds_[dim].lo; + result = result - floor(temp / box_width + 1) * box_width; + } + } else { + result = temp; + } + + return result; +} + +template +double DHrectBound::MinDelta(const DHrectBound& other, double box_width, int dim) const { + double result = -0.5 * box_width; + double temp = other.bounds_[dim].hi - bounds_[dim].lo; + temp = temp - floor(temp / box_width) * box_width; + if (temp > box_width / 2) { + temp = other.bounds_[dim].hi - bounds_[dim].hi; + temp -= floor(temp / box_width) * box_width; + if (temp > box_width / 2) + result = temp - box_width; + + } else { + temp = other.bounds_[dim].hi - bounds_[dim].hi; + result = temp - floor(temp / box_width) * box_width; + } + + return result; +} + +/** + * Calculates minimum and maximum bound-to-bound squared distance. + */ +template +DRange DHrectBound::RangeDistanceSq(const DHrectBound& other) const { + double sum_lo = 0; + double sum_hi = 0; + + DEBUG_SAME_SIZE(dim_, other.dim_); + + double v1, v2, v_lo, v_hi; + for (index_t d = 0; d < dim_; d++) { + v1 = other.bounds_[d].lo - bounds_[d].hi; + v2 = bounds_[d].lo - other.bounds_[d].hi; + // one of v1 or v2 is negative + if(v1 >= v2) { + v_hi = -v2; // make it nonnegative + v_lo = (v1 > 0) ? v1 : 0; // force to be 0 if negative + } else { + v_hi = -v1; // make it nonnegative + v_lo = (v2 > 0) ? v2 : 0; // force to be 0 if negative + } + + sum_lo += pow(v_lo, (double) t_pow); + sum_hi += pow(v_hi, (double) t_pow); + } + + return DRange(pow(sum_lo, 2.0 / (double) t_pow), + pow(sum_hi, 2.0 / (double) t_pow)); +} + +/** + * Calculates minimum and maximum bound-to-point squared distance. + */ +template +DRange DHrectBound::RangeDistanceSq(const vec& point) const { + double sum_lo = 0; + double sum_hi = 0; + + DEBUG_SAME_SIZE(point.n_elem, dim_); + + double v1, v2, v_lo, v_hi; + for(index_t d = 0; d < dim_; d++) { + v1 = bounds_[d].lo - point[d]; + v2 = point[d] - bounds_[d].hi; + // one of v1 or v2 is negative + if(v1 >= 0) { + v_hi = -v2; + v_lo = v1; + } else { + v_hi = -v1; + v_lo = v2; + } + + sum_lo += pow(v_lo, (double) t_pow); + sum_hi += pow(v_hi, (double) t_pow); + } + + return DRange(pow(sum_lo, 2.0 / (double) t_pow), + pow(sum_hi, 2.0 / (double) t_pow)); +} + +/** + * Calculates closest-to-their-midpoint bounding box distance, + * i.e. calculates their midpoint and finds the minimum box-to-point + * distance. + * + * Equivalent to: + * + * other.CalcMidpoint(&other_midpoint) + * return MinDistanceSqToPoint(other_midpoint) + * + */ +template +double DHrectBound::MinToMidSq(const DHrectBound& other) const { + double sum = 0; + + DEBUG_SAME_SIZE(dim_, other.dim_); + + for (index_t d = 0; d < dim_; d++) { + double v = other.bounds_[d].mid(); + double v1 = bounds_[d].lo - v; + double v2 = v - bounds_[d].hi; + + v = (v1 + fabs(v1)) + (v2 + fabs(v2)); + + sum += pow(v, (double) t_pow); // v is non-negative + } + + return pow(sum, 2.0 / (double) t_pow) / 4.0; +} + +/** + * Computes minimax distance, where the other node is trying to avoid me. + */ +template +double DHrectBound::MinimaxDistanceSq(const DHrectBound& other) const { + double sum = 0; + + DEBUG_SAME_SIZE(dim_, other.dim_); + + for(index_t d = 0; d < dim_; d++) { + double v1 = other.bounds_[d].hi - bounds_[d].hi; + double v2 = bounds_[d].lo - other.bounds_[d].lo; + double v = std::max(v1, v2); + v = (v + fabs(v)); /* truncate negatives to zero */ + sum += pow(v, (double) t_pow); // v is non-negative + } + + return pow(sum, 2.0 / (double) t_pow) / 4.0; +} + +/** + * Calculates midpoint-to-midpoint bounding box distance. + */ +template +double DHrectBound::MidDistanceSq(const DHrectBound& other) const { + double sum = 0; + const DRange *a = this->bounds_; + const DRange *b = other.bounds_; + + DEBUG_SAME_SIZE(dim_, other.dim_); + + for (index_t d = 0; d < dim_; d++) { + // take the midpoint of each dimension (left multiplied by two for + // calculation speed) and subtract from each other, then raise to t_pow + sum += pow(fabs(bounds_[d].hi + bounds_[d].lo - other.bounds_[d].hi - other.bounds_[d].lo), (double) t_pow); + } + + // take t_pow/2'th root and divide by constant of 4 (leftover from previous + // step) + return pow(sum, 2.0 / (double) t_pow) / 4.0; +} + +/** + * Expands this region to include a new point. + */ +template +DHrectBound& DHrectBound::operator|=(const vec& vector) { + DEBUG_SAME_SIZE(vector.n_elem, dim_); + + for (index_t i = 0; i < dim_; i++) { + bounds_[i] |= vector[i]; + } + + return *this; +} + +/** + * Expands this region to encompass another bound. + */ +template +DHrectBound& DHrectBound::operator|=(const DHrectBound& other) { + DEBUG_SAME_SIZE(other.dim_, dim_); + + for (index_t i = 0; i < dim_; i++) { + bounds_[i] |= other.bounds_[i]; + } + + return *this; +} + + +/** + * Expand this bounding box to encompass another point. Done to + * minimize added volume in periodic coordinates. + */ +template +DHrectBound& DHrectBound::Add(const vec& other, const vec& size) { + DEBUG_SAME_SIZE(other.n_elem, dim_); + // Catch case of uninitialized bounds + if (bounds_[0].hi < 0){ + for (index_t i = 0; i < dim_; i++){ + bounds_[i] |= other[i]; + } + } + + for (index_t i= 0; i < dim_; i++){ + double ah, al; + ah = bounds_[i].hi - other[i]; + al = bounds_[i].lo - other[i]; + ah = ah - floor(ah / size[i]) * size[i]; + al = al - floor(al / size[i]) * size[i]; + if (ah < al) { + if (size[i] - ah < al) { + bounds_[i].hi = other[i]; + } else { + bounds_[i].lo = other[i]; + } + } + } + return *this; +} + +/** + * Expand this bounding box in periodic coordinates, minimizing added volume. + */ +template +DHrectBound& DHrectBound::Add(const DHrectBound& other, const vec& size){ + if (bounds_[0].hi < 0){ + for (index_t i = 0; i < dim_; i++){ + bounds_[i] |= other.bounds_[i]; + } + } + + for (index_t i = 0; i < dim_; i++) { + double ah, al, bh, bl; + ah = bounds_[i].hi; + al = bounds_[i].lo; + bh = other.bounds_[i].hi; + bl = other.bounds_[i].lo; + ah = ah - al; + bh = bh - al; + bl = bl - al; + ah = ah - floor(ah / size[i]) * size[i]; + bh = bh - floor(bh / size[i]) * size[i]; + bl = bl - floor(bl / size[i]) * size[i]; + + if (((bh > ah) & (bh < bl | ah > bl )) || + (bh >= bl & bl > ah & bh < ah -bl + size[i])){ + bounds_[i].hi = other.bounds_[i].hi; + } + + if (bl > ah && ((bl > bh) || (bh >= ah -bl + size[i]))){ + bounds_[i].lo = other.bounds_[i].lo; + } + + if (unlikely(ah > bl & bl > bh)){ + bounds_[i].lo = 0; + bounds_[i].hi = size[i]; + } + } + + return *this; +} + +#endif diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/kdtree.h b/fastlib/branches/fastlib-stl/fastlib/tree/kdtree.h index 39dc6515e0..ca38ea16a6 100644 --- a/fastlib/branches/fastlib-stl/fastlib/tree/kdtree.h +++ b/fastlib/branches/fastlib-stl/fastlib/tree/kdtree.h @@ -47,15 +47,14 @@ #include "spacetree.h" #include "bounds.h" -//#include "spacetree.h" -//#include "bounds.h" #include #include "../fx/fx.h" #include "kdtree_impl.h" -//#include "kdtree_impl.h" + +#include /** * Regular pointer-style trees (as opposed to THOR trees). @@ -66,127 +65,162 @@ namespace tree { * * @experimental * - * This requires you to pass in two unitialized ArrayLists which will contain - * index mappings so you can account for the re-ordering of the matrix. - * (By unitialized I mean don't call Init on it) + * This requires you to pass in two vectors which will contain index mappings + * so you can account for the re-ordering of the matrix. * * @param matrix data where each column is a point, WHICH WILL BE RE-ORDERED + * @param split_dimensions ordering of the dimensions that we should split on; + * the first element in this vector will be the first element we split + * on, and so on * @param leaf_size the maximum points in a leaf - * @param old_from_new pointer to an unitialized arraylist; it will map - * new indices to original - * @param new_from_old pointer to an unitialized arraylist; it will map - * original indexes to new indices + * @param old_from_new vector that will map new indices to original + * @param new_from_old vector that will map original indexes to new indices */ - - template - TKdTree *MakeKdTreeMidpointSelective(GenMatrix& matrix, - const Vector& split_dimensions, - index_t leaf_size, - std::vector *old_from_new = NULL, - std::vector *new_from_old = NULL) { + TKdTree *MakeKdTreeMidpointSelective(arma::Mat& matrix, + const arma::uvec& split_dimensions, + index_t leaf_size, + arma::Col& old_from_new, + arma::Col& new_from_old) { TKdTree *node = new TKdTree(); - index_t *old_from_new_ptr; - if (old_from_new) { - old_from_new->reserve(matrix.n_cols()); + old_from_new.set_size(matrix.n_cols); + for (index_t i = 0; i < matrix.n_cols; i++) + old_from_new[i] = i; - for (index_t i = 0; i < matrix.n_cols(); i++) { - (*old_from_new)[i] = i; - } - - old_from_new_ptr = old_from_new->begin(); - } else { - old_from_new_ptr = NULL; - } - - node->Init(0, matrix.n_cols()); - node->bound().Init(split_dimensions.length()); + node->Init(0, matrix.n_cols); + node->bound().SetSize(split_dimensions.n_elem); tree_kdtree_private::SelectFindBoundFromMatrix(matrix, split_dimensions, - 0, matrix.n_cols(), &node->bound()); + 0, matrix.n_cols, &node->bound()); tree_kdtree_private::SelectSplitKdTreeMidpoint(matrix, split_dimensions, - node, leaf_size, old_from_new_ptr); + node, leaf_size, old_from_new); - if (new_from_old) { - new_from_old->reserve(matrix.n_cols()); - for (index_t i = 0; i < matrix.n_cols(); i++) { - (*new_from_old)[(*old_from_new)[i]] = i; + if(new_from_old) { + new_from_old.set_size(matrix.n_cols); + for (index_t i = 0; i < matrix.n_cols; i++) { + new_from_old[old_from_new[i]] = i; } - } + } + return node; + } + + /** + * Creates a KD tree from data, splitting on the midpoint. + * + * @experimental + * + * This requires you to pass in one vector that will map new indices to the + * original. It does not bother keeping track of a map of the original + * indices to the new indices. + * + * @param matrix data where each column is a point, WHICH WILL BE RE-ORDERED + * @param split_dimensions ordering of the dimensions that we should split on; + * the first element in this vector will be the first element we split + * on, and so on + * @param leaf_size the maximum points in a leaf + * @param old_from_new vector that will map new indices to original + */ + template + TKdTree *MakeKdTreeMidpointSelective(arma::Mat& matrix, + const arma::uvec& split_dimensions, + index_t leaf_size, + arma::Col& old_from_new) { + TKdTree *node = new TKdTree(); + + old_from_new.set_size(matrix.n_cols); + for (index_t i = 0; i < matrix.n_cols; i++) + old_from_new[i] = i; + + node->Init(0, matrix.n_cols); + node->bound().SetSize(split_dimensions.n_elem); + tree_kdtree_private::SelectFindBoundFromMatrix(matrix, split_dimensions, + 0, matrix.n_cols, &node->bound()); + + tree_kdtree_private::SelectSplitKdTreeMidpoint(matrix, split_dimensions, + node, leaf_size, old_from_new); + + return node; + } + + /** + * Creates a KD tree from data, splitting on the midpoint. + * + * @experimental + * + * This does not keep track of a map of the old indices to the new indices. + * + * @param matrix data where each column is a point, WHICH WILL BE RE-ORDERED + * @param split_dimensions ordering of the dimensions that we should split on; + * the first element in this vector will be the first element we split + * on, and so on + * @param leaf_size the maximum points in a leaf + */ + template + TKdTree *MakeKdTreeMidpointSelective(arma::Mat& matrix, + const arma::uvec& split_dimensions, + index_t leaf_size) { + TKdTree *node = new TKdTree(); + + node->Init(0, matrix.n_cols); + node->bound().SetSize(split_dimensions.n_elem); + tree_kdtree_private::SelectFindBoundFromMatrix(matrix, split_dimensions, + 0, matrix.n_cols, &node->bound()); + + tree_kdtree_private::SelectSplitKdTreeMidpoint(matrix, split_dimensions, + node, leaf_size); + return node; } template - TKdTree *MakeKdTreeMidpointSelective(GenMatrix& matrix, - const Vector& split_dimensions, - index_t leaf_size, - GenVector *old_from_new = NULL, - GenVector *new_from_old = NULL) { - TKdTree *node = new TKdTree(); - index_t *old_from_new_ptr; - - if (old_from_new) { - old_from_new->Init(matrix.n_cols()); - - for (index_t i = 0; i < matrix.n_cols(); i++) { - (*old_from_new)[i] = i; - } - - old_from_new_ptr = old_from_new->ptr(); - } else { - old_from_new_ptr = NULL; - } - - node->Init(0, matrix.n_cols()); - node->bound().Init(split_dimensions.length()); - tree_kdtree_private::SelectFindBoundFromMatrix(matrix, split_dimensions, - 0, matrix.n_cols(), &node->bound()); - - tree_kdtree_private::SelectSplitKdTreeMidpoint(matrix, split_dimensions, - node, leaf_size, old_from_new_ptr); - - if (new_from_old) { - new_from_old->Init(matrix.n_cols()); - for (index_t i = 0; i < matrix.n_cols(); i++) { - (*new_from_old)[(*old_from_new)[i]] = i; - } - } - return node; - } - - - template - TKdTree *MakeKdTreeMidpoint(GenMatrix& matrix, index_t leaf_size, - std::vector *old_from_new = NULL, - std::vector *new_from_old = NULL) { - Vector split_dimensions; - split_dimensions.Init(matrix.n_rows()); - int i; - for (i = 0; i < matrix.n_rows(); i++){ + TKdTree *MakeKdTreeMidpoint(arma::Mat& matrix, + index_t leaf_size, + arma::Col& old_from_new, + arma::Col& new_from_old) { + // create vector of dimensions that we will split on + // by default we'll just split the first dimension first, and so on + arma::uvec split_dimensions(matrix.n_rows); + for(int i = 0; i < matrix.n_rows; i++) split_dimensions[i] = i; - } + TKdTree *result; - result = MakeKdTreeMidpointSelective(matrix, split_dimensions, - leaf_size, old_from_new, new_from_old); + result = MakeKdTreeMidpointSelective(matrix, + split_dimensions, leaf_size, old_from_new, new_from_old); + return result; } template - TKdTree *MakeKdTreeMidpoint(GenMatrix& matrix, - index_t leaf_size, - GenVector *old_from_new = NULL, - GenVector *new_from_old = NULL) { - Vector split_dimensions; - split_dimensions.Init(matrix.n_rows()); - int i; - for (i = 0; i < matrix.n_rows(); i++){ + TKdTree *MakeKdTreeMidpoint(arma::Mat& matrix, + index_t leaf_size, + arma::Col& old_from_new) { + // create vector of dimensions that we will split on + // by default we'll just split the first dimension first, and so on + arma::uvec split_dimensions(matrix.n_rows); + for(int i = 0; i < matrix.n_rows; i++) split_dimensions[i] = i; - } + TKdTree *result; result = MakeKdTreeMidpointSelective(matrix, - split_dimensions, - leaf_size, old_from_new, new_from_old); + split_dimensions, leaf_size, old_from_new); + + return result; + } + + template + TKdTree *MakeKdTreeMidpoint(arma::Mat& matrix, + index_t leaf_size) { + // create vector of dimensions that we will split on + // by default we'll just split the first dimension first, and so on + arma::uvec split_dimensions(matrix.n_rows); + for(int i = 0; i < matrix.n_rows; i++) + split_dimensions[i] = i; + + TKdTree *result; + result = MakeKdTreeMidpointSelective(matrix, + split_dimensions, leaf_size); + return result; } @@ -208,7 +242,7 @@ namespace tree { * @code * MyKdTree *q_tree; * Matrix q_matrix; - * std::vector q_permutation; + * arma::Col q_permutation; * LoadKdTree(fx_submodule(NULL, "q", "q"), &q_matrix, &q_tree, * &q_permutation); * @endcode @@ -229,9 +263,10 @@ namespace tree { * @return SUCCESS_PASS or SUCCESS_FAIL */ template - success_t LoadKdTree(datanode *module, - GenMatrix *matrix, TKdTree **tree_pp, - std::vector *old_from_new) { + success_t LoadKdTree(datanode* module, + arma::Mat& matrix, + TKdTree** tree_pp, + arma::Col& old_from_new) { const char *type = fx_param_str(module, "type", "text"); const char *fname = fx_param_str(module, "", NULL); success_t success = SUCCESS_PASS; @@ -248,7 +283,7 @@ namespace tree { fx_timer_start(module, "make_tree"); *tree_pp = MakeKdTreeMidpoint( - *matrix, leaflen, old_from_new); + matrix, leaflen, old_from_new); fx_timer_stop(module, "make_tree"); } fx_timer_stop(module, "load"); diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/kdtree_impl.cc b/fastlib/branches/fastlib-stl/fastlib/tree/kdtree_impl.cc new file mode 100644 index 0000000000..1b8162dda6 --- /dev/null +++ b/fastlib/branches/fastlib-stl/fastlib/tree/kdtree_impl.cc @@ -0,0 +1,21 @@ +/** + * @file kdtree_impl.cc + * + * Specialized functions. + */ +#include "../fastlib.h" +#include "../la/matrix.h" +#include "kdtree.h" + +#include + +/*** + * Specialized MakeBoundVector function written for arma::vec instead of + * GenVector. + */ +void tree_kdtree_private::MakeBoundVector(const arma::vec& point, + const arma::uvec& bound_dimensions, + arma::vec& bound_vector) { + for(int i = 0; i < bound_dimensions.n_elem; i++) + bound_vector[i] = point[(int) bound_dimensions[i]]; +} diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/kdtree_impl.h b/fastlib/branches/fastlib-stl/fastlib/tree/kdtree_impl.h index 1be27e531c..26ba55a72f 100644 --- a/fastlib/branches/fastlib-stl/fastlib/tree/kdtree_impl.h +++ b/fastlib/branches/fastlib-stl/fastlib/tree/kdtree_impl.h @@ -31,70 +31,92 @@ */ /* Implementation for the regular pointer-style kd-tree builder. */ +#ifndef TREE_KDTREE_IMPL_H +#define TREE_KDTREE_IMPL_H + +#include "../base/arma_compat.h" + namespace tree_kdtree_private { - template - void MakeBoundVector(const GenVector& point, - const Vector& bound_dimensions, GenVector* bound_vector){ - int i; - for (i = 0; i < bound_dimensions.length(); i++){ - (*bound_vector)[i] = point[(int)bound_dimensions[i]]; - } - } - + void MakeBoundVector(const arma::vec& point, + const arma::uvec& bound_dimensions, + arma::vec& bound_vector); template - void SelectFindBoundFromMatrix(const GenMatrix& matrix, - const Vector& split_dimensions, index_t first, index_t count, - TBound *bounds){ + void SelectFindBoundFromMatrix(const arma::Mat& matrix, + const arma::uvec& split_dimensions, + index_t first, + index_t count, + TBound *bounds) { index_t end = first + count; - for (index_t i = first; i < end; i++) { - GenVector col; - matrix.MakeColumnVector(i, &col); - if (split_dimensions.length() == matrix.n_rows()){ - *bounds |= col; + for(index_t i = first; i < end; i++) { + if (split_dimensions.n_elem == matrix.n_rows){ + *bounds |= matrix.col(i); } else { - GenVector sub_col; - sub_col.Init(split_dimensions.length()); - MakeBoundVector(col, split_dimensions, &sub_col); + arma::vec sub_col(split_dimensions.n_elem); + MakeBoundVector(matrix.col(i), split_dimensions, sub_col); *bounds |= sub_col; } } } template - void FindBoundFromMatrix(const GenMatrix& matrix, - index_t first, index_t count, TBound *bounds){ - Vector split_dimensions; - split_dimensions.Init(matrix.n_rows()); - int i; - for (i = 0; i < matrix.n_rows(); i++){ + void FindBoundFromMatrix(const arma::Mat& matrix, + index_t first, + index_t count, + TBound *bounds) { + arma::uvec split_dimensions(matrix.n_rows); + for(int i = 0; i < matrix.n_rows; i++) split_dimensions[i] = i; - } + SelectFindBoundFromMatrix(matrix, split_dimensions, first, count, bounds); } - template - index_t MatrixPartition( - Matrix& matrix, index_t dim, double splitvalue, - index_t first, index_t count, - TBound* left_bound, TBound* right_bound, - index_t *old_from_new) { - Vector split_dimensions; - split_dimensions.Init(matrix.n_rows()); - int i; - for (i = 0; i < matrix.n_rows(); i++){ - split_dimensions[i] = i; - } - index_t split_point = SelectMatrixPartition(matrix, split_dimensions, - dim, splitvalue, first, count, left_bound, right_bound, old_from_new); - return split_point; - } + template + index_t MatrixPartition(arma::mat& matrix, + index_t dim, + double splitvalue, + index_t first, + index_t count, + TBound& left_bound, + TBound& right_bound, + arma::Col& old_from_new) { + // we will split dimensions in a very simple order: first dim first + arma::uvec split_dimensions(matrix.n_rows); + for(int i = 0; i < matrix.n_rows; i++) + split_dimensions[i] = i; + index_t split_point = SelectMatrixPartition(matrix, split_dimensions, + dim, splitvalue, first, count, left_bound, right_bound, old_from_new); + return split_point; + } + + template + index_t MatrixPartition(arma::mat& matrix, + index_t dim, + double splitvalue, + index_t first, + index_t count, + TBound& left_bound, + TBound& right_bound) { + // we will split dimensions in a very simple order: first dim first + arma::uvec split_dimensions(matrix.n_rows); + for(int i = 0; i < matrix.n_rows; i++) + split_dimensions[i] = i; + + index_t split_point = SelectMatrixPartition(matrix, split_dimensions, + dim, splitvalue, first, count, left_bound, right_bound); + return split_point; + } + template - index_t SelectMatrixPartition(GenMatrix& matrix, - const Vector& split_dimensions, index_t dim, double splitvalue, - index_t first, index_t count, TBound* left_bound, TBound* right_bound, - index_t *old_from_new) { + index_t SelectMatrixPartition(arma::Mat& matrix, + const arma::uvec& split_dimensions, + index_t dim, + double splitvalue, + index_t first, + index_t count, + TBound& left_bound, + TBound& right_bound) { index_t left = first; index_t right = first + count - 1; @@ -105,30 +127,26 @@ namespace tree_kdtree_private { * everything > right is correct */ for (;;) { - while (matrix.get(dim, left) < splitvalue && likely(left <= right)) { - GenVector left_vector; - matrix.MakeColumnVector(left, &left_vector); - if (split_dimensions.length() == matrix.n_rows()){ - *left_bound |= left_vector; + while (matrix(dim, left) < splitvalue && likely(left <= right)) { + arma::vec left_vector = matrix.col(left); + if (split_dimensions.n_elem == matrix.n_rows) { + left_bound |= left_vector; } else { - GenVector sub_left_vector; - sub_left_vector.Init(split_dimensions.length()); - MakeBoundVector(left_vector, split_dimensions, &sub_left_vector); - *left_bound |= sub_left_vector; + arma::vec sub_left_vector(split_dimensions.n_elem); + MakeBoundVector(left_vector, split_dimensions, sub_left_vector); + left_bound |= sub_left_vector; } left++; } - while (matrix.get(dim, right) >= splitvalue && likely(left <= right)) { - GenVector right_vector; - matrix.MakeColumnVector(right, &right_vector); - if (split_dimensions.length() == matrix.n_rows()){ - *right_bound |= right_vector; + while (matrix(dim, right) >= splitvalue && likely(left <= right)) { + arma::vec right_vector = matrix.col(right); + if (split_dimensions.n_elem == matrix.n_rows) { + right_bound |= right_vector; } else { - GenVector sub_right_vector; - sub_right_vector.Init(split_dimensions.length()); - MakeBoundVector(right_vector, split_dimensions, &sub_right_vector); - *right_bound |= sub_right_vector; + arma::vec sub_right_vector(split_dimensions.n_elem); + MakeBoundVector(right_vector, split_dimensions, sub_right_vector); + right_bound |= sub_right_vector; } right--; } @@ -138,46 +156,115 @@ namespace tree_kdtree_private { break; } - GenVector left_vector; - GenVector right_vector; + // swap left and right vector + matrix.swap_cols(left, right); - matrix.MakeColumnVector(left, &left_vector); - matrix.MakeColumnVector(right, &right_vector); + arma::vec left_vector = matrix.col(left); + arma::vec right_vector = matrix.col(right); - left_vector.SwapValues(&right_vector); - - if (split_dimensions.length() == matrix.n_rows()){ - *left_bound |= left_vector; + if (split_dimensions.n_elem == matrix.n_rows) { + left_bound |= left_vector; } else { - GenVector sub_left_vector; - sub_left_vector.Init(split_dimensions.length()); - MakeBoundVector(left_vector, split_dimensions, &sub_left_vector); - *left_bound |= sub_left_vector; + arma::vec sub_left_vector(split_dimensions.n_elem); + MakeBoundVector(left_vector, split_dimensions, sub_left_vector); + left_bound |= sub_left_vector; } - if (split_dimensions.length() == matrix.n_rows()){ - *right_bound |= right_vector; + if (split_dimensions.n_elem == matrix.n_rows){ + right_bound |= right_vector; } else { - GenVector sub_right_vector; - sub_right_vector.Init(split_dimensions.length()); - MakeBoundVector(right_vector, split_dimensions, &sub_right_vector); - *right_bound |= sub_right_vector; + arma::vec sub_right_vector(split_dimensions.n_elem); + MakeBoundVector(right_vector, split_dimensions, sub_right_vector); + right_bound |= sub_right_vector; } - - if (old_from_new) { - index_t t = old_from_new[left]; - old_from_new[left] = old_from_new[right]; - old_from_new[right] = t; + DEBUG_ASSERT(left <= right); + right--; + } + + DEBUG_ASSERT(left == right + 1); + + return left; + } + + template + index_t SelectMatrixPartition(arma::Mat& matrix, + const arma::uvec& split_dimensions, + index_t dim, + double splitvalue, + index_t first, + index_t count, + TBound& left_bound, + TBound& right_bound, + arma::Col& old_from_new) { + + index_t left = first; + index_t right = first + count - 1; + + /* At any point: + * + * everything < left is correct + * everything > right is correct + */ + for (;;) { + while (matrix(dim, left) < splitvalue && likely(left <= right)) { + arma::vec left_vector = matrix.col(left); + if (split_dimensions.n_elem == matrix.n_rows) { + left_bound |= left_vector; + } else { + arma::vec sub_left_vector(split_dimensions.n_elem); + MakeBoundVector(left_vector, split_dimensions, sub_left_vector); + left_bound |= sub_left_vector; + } + left++; } + while (matrix(dim, right) >= splitvalue && likely(left <= right)) { + arma::vec right_vector = matrix.col(right); + if (split_dimensions.n_elem == matrix.n_rows) { + right_bound |= right_vector; + } else { + arma::vec sub_right_vector(split_dimensions.n_elem); + MakeBoundVector(right_vector, split_dimensions, sub_right_vector); + right_bound |= sub_right_vector; + } + right--; + } + + if (unlikely(left > right)) { + /* left == right + 1 */ + break; + } + + // swap left and right vector + matrix.swap_cols(left, right); + + arma::vec left_vector = matrix.col(left); + arma::vec right_vector = matrix.col(right); + + if (split_dimensions.n_elem == matrix.n_rows) { + left_bound |= left_vector; + } else { + arma::vec sub_left_vector(split_dimensions.n_elem); + MakeBoundVector(left_vector, split_dimensions, sub_left_vector); + left_bound |= sub_left_vector; + } + + if (split_dimensions.n_elem == matrix.n_rows){ + right_bound |= right_vector; + } else { + arma::vec sub_right_vector(split_dimensions.n_elem); + MakeBoundVector(right_vector, split_dimensions, sub_right_vector); + right_bound |= sub_right_vector; + } + + // update indices + index_t t = old_from_new[left]; + old_from_new[left] = old_from_new[right]; + old_from_new[right] = t; + DEBUG_ASSERT(left <= right); right--; - - // this conditional is always true, I belueve - //if (likely(left <= right)) { - // right--; - //} } DEBUG_ASSERT(left == right + 1); @@ -186,24 +273,38 @@ namespace tree_kdtree_private { } template - void SplitKdTreeMidpoint(GenMatrix& matrix, - TKdTree *node, index_t leaf_size, index_t *old_from_new){ + void SplitKdTreeMidpoint(arma::Mat& matrix, + TKdTree *node, + index_t leaf_size, + arma::Col& old_from_new) { - Vector split_dimensions; - split_dimensions.Init(matrix.n_rows()); - int i; - for (i = 0; i < matrix.n_rows(); i++){ + arma::uvec split_dimensions(matrix.n_rows); + for(int i = 0; i < matrix.n_rows; i++) split_dimensions[i] = i; - } + SelectSplitKdTreeMidpoint(matrix, split_dimensions, node, leaf_size, old_from_new); } + template + void SplitKdTreeMidpoint(arma::Mat& matrix, + TKdTree *node, + index_t leaf_size) { + + arma::uvec split_dimensions(matrix.n_rows); + for(int i = 0; i < matrix.n_rows; i++) + split_dimensions[i] = i; + + SelectSplitKdTreeMidpoint(matrix, split_dimensions, node, + leaf_size); + } template - void SelectSplitKdTreeMidpoint(GenMatrix& matrix, - const Vector& split_dimensions, TKdTree *node, index_t leaf_size, - index_t *old_from_new) { + void SelectSplitKdTreeMidpoint(arma::Mat& matrix, + const arma::uvec& split_dimensions, + TKdTree *node, + index_t leaf_size, + arma::Col& old_from_new) { TKdTree *left = NULL; TKdTree *right = NULL; @@ -211,43 +312,43 @@ namespace tree_kdtree_private { SelectFindBoundFromMatrix(matrix, split_dimensions, node->begin(), node->count(), &node->bound()); - if (node->count() > leaf_size) { + if(node->count() > leaf_size) { index_t split_dim = BIG_BAD_NUMBER; double max_width = -1; - for (index_t d = 0; d < split_dimensions.length(); d++) { - double w = node->bound().get(d).width(); + for (index_t d = 0; d < split_dimensions.n_elem; d++) { + double w = node->bound()[d].width(); - if (w > max_width) { + if (w > max_width) { max_width = w; split_dim = d; } } - double split_val = node->bound().get(split_dim).mid(); + double split_val = node->bound()[split_dim].mid(); if (max_width == 0) { // Okay, we can't do any splitting, because all these points are the // same. We have to give up. } else { left = new TKdTree(); - left->bound().Init(split_dimensions.length()); + left->bound().SetSize(split_dimensions.n_elem); right = new TKdTree(); - right->bound().Init(split_dimensions.length()); + right->bound().SetSize(split_dimensions.n_elem); index_t split_col = SelectMatrixPartition(matrix, split_dimensions, - (int)split_dimensions[split_dim], split_val, + (int) split_dimensions[split_dim], split_val, node->begin(), node->count(), - &left->bound(), &right->bound(), + left->bound(), right->bound(), old_from_new); VERBOSE_MSG(3.0,"split (%d,[%d],%d) dim %d on %f (between %f, %f)", node->begin(), split_col, - node->begin() + node->count(), (int)split_dimensions[split_dim], + node->begin() + node->count(), (int) split_dimensions[split_dim], split_val, - node->bound().get(split_dim).lo, - node->bound().get(split_dim).hi); + node->bound()[split_dim].lo, + node->bound()[split_dim].hi); left->Init(node->begin(), split_col - node->begin()); right->Init(split_col, node->begin() + node->count() - split_col); @@ -264,4 +365,69 @@ namespace tree_kdtree_private { node->set_children(matrix, left, right); } + + template + void SelectSplitKdTreeMidpoint(arma::Mat& matrix, + const arma::uvec& split_dimensions, + TKdTree *node, + index_t leaf_size) { + + TKdTree *left = NULL; + TKdTree *right = NULL; + + SelectFindBoundFromMatrix(matrix, split_dimensions, node->begin(), + node->count(), &node->bound()); + + if(node->count() > leaf_size) { + index_t split_dim = BIG_BAD_NUMBER; + double max_width = -1; + + for (index_t d = 0; d < split_dimensions.n_elem; d++) { + double w = node->bound()[d].width(); + + if (w > max_width) { + max_width = w; + split_dim = d; + } + } + + double split_val = node->bound()[split_dim].mid(); + + if (max_width == 0) { + // Okay, we can't do any splitting, because all these points are the + // same. We have to give up. + } else { + left = new TKdTree(); + left->bound().SetSize(split_dimensions.n_elem); + + right = new TKdTree(); + right->bound().SetSize(split_dimensions.n_elem); + + index_t split_col = SelectMatrixPartition(matrix, split_dimensions, + (int) split_dimensions[split_dim], split_val, + node->begin(), node->count(), + left->bound(), right->bound()); + + VERBOSE_MSG(3.0,"split (%d,[%d],%d) dim %d on %f (between %f, %f)", + node->begin(), split_col, + node->begin() + node->count(), (int) split_dimensions[split_dim], + split_val, + node->bound()[split_dim].lo, + node->bound()[split_dim].hi); + + left->Init(node->begin(), split_col - node->begin()); + right->Init(split_col, node->begin() + node->count() - split_col); + + // This should never happen if max_width > 0 + DEBUG_ASSERT(left->count() != 0 && right->count() != 0); + + SelectSplitKdTreeMidpoint(matrix, split_dimensions, left, leaf_size); + SelectSplitKdTreeMidpoint(matrix, split_dimensions, right, leaf_size); + } + } + + node->set_children(matrix, left, right); + } }; + +#endif diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/lmetric.h b/fastlib/branches/fastlib-stl/fastlib/tree/lmetric.h new file mode 100644 index 0000000000..4794ff87e1 --- /dev/null +++ b/fastlib/branches/fastlib-stl/fastlib/tree/lmetric.h @@ -0,0 +1,83 @@ +/* MLPACK 0.2 + * + * Copyright (c) 2008, 2009 Alexander Gray, + * Garry Boyer, + * Ryan Riegel, + * Nikolaos Vasiloglou, + * Dongryeol Lee, + * Chip Mappus, + * Nishant Mehta, + * Hua Ouyang, + * Parikshit Ram, + * Long Tran, + * Wee Chin Wong + * + * Copyright (c) 2008, 2009 Georgia Institute of Technology + * + * This program 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 2 of the + * License, or (at your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA + * 02110-1301, USA. + */ +/** + * @file tree/lmetric.h + * + * Bounds that are useful for binary space partitioning trees. + * + * This file describes the LMetric policy class. + * + * @experimental + */ + +#ifndef TREE_LMETRIC_H +#define TREE_LMETRIC_H + +#include "../la/la.h" +#include "../la/matrix.h" + +#include +#include + +/** + * An L_p metric for vector spaces. + * + * A generic Metric class should simply compute the distance between + * two points. An LMetric operates for integer powers on vector spaces. + */ +template +class LMetric { + public: + /** + * Computes the distance metric between two points. + */ + static double Distance(const arma::vec& a, const arma::vec& b) { + // subtract b from a elementwise, and then raise to t_pow; + // then sum all that and take the t_pow'th root of it + return pow(accu(pow((a - b), t_pow)), 1.0 / (double) t_pow); + } + + /** + * Computes the distance metric between two points, raised to a + * particular power. + * + * This might be faster so that you could get, for instance, squared + * L2 distance. + */ + template + static double PowDistance(const arma::vec& a, const arma::vec& b) { + // accu() sums all elements of a vector + return pow(accu(pow((a - b), t_pow)), (double) t_result_pow / (double) t_pow); + } +}; + +#endif diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/spacetree.h b/fastlib/branches/fastlib-stl/fastlib/tree/spacetree.h index c72262ea4f..69767bc7f7 100644 --- a/fastlib/branches/fastlib-stl/fastlib/tree/spacetree.h +++ b/fastlib/branches/fastlib-stl/fastlib/tree/spacetree.h @@ -42,7 +42,9 @@ #include "../base/base.h" #include "statistic.h" -//#include "statistic.h" + +#include +#include "../base/arma_compat.h" /** * A binary space partitioning tree, such as KD or ball tree. @@ -167,7 +169,7 @@ class BinarySpaceTree { /** * Used only when constructing the tree. */ - void set_children(const Dataset& data, + void set_children(const TDataset& data, BinarySpaceTree *left_in, BinarySpaceTree *right_in) { left_ = left_in; right_ = right_in; @@ -181,6 +183,27 @@ class BinarySpaceTree { } } + /** + * Function specialization for arma transition. + * Used only when constructing the tree. + * This assumes TDataset is Dataset. + * TODO: remove this function + */ + void set_children(const arma::mat& data, + BinarySpaceTree *left_in, BinarySpaceTree *right_in) { + left_ = left_in; + right_ = right_in; + + // make a fake matrix + Matrix tmp; + arma_compat::armaToMatrix(data, tmp); + if(!is_leaf()) + stat_.Init(tmp, begin_, count_, left_->stat_, right_->stat_); + else + stat_.Init(tmp, begin_, count_); + + } + const Bound& bound() const { return bound_; } diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/statistic.h b/fastlib/branches/fastlib-stl/fastlib/tree/statistic.h index ab04b1f68e..aa28a5e211 100644 --- a/fastlib/branches/fastlib-stl/fastlib/tree/statistic.h +++ b/fastlib/branches/fastlib-stl/fastlib/tree/statistic.h @@ -50,24 +50,22 @@ */ template class EmptyStatistic { - public: - EmptyStatistic() {} - ~EmptyStatistic() {} + public: + EmptyStatistic() {} + ~EmptyStatistic() {} - /** - * Initializes by taking statistics on raw data. - */ - void Init(const TDataset& dataset, index_t start, index_t count) { - } + /** + * Initializes by taking statistics on raw data. + */ + void Init(const TDataset& dataset, index_t start, index_t count) { } - /** - * Initializes by combining statistics of two partitions. - * - * This lets you build fast bottom-up statistics when building trees. - */ - void Init(const TDataset& dataset, index_t start, index_t count, - const EmptyStatistic& left_stat, const EmptyStatistic& right_stat) { - } + /** + * Initializes by combining statistics of two partitions. + * + * This lets you build fast bottom-up statistics when building trees. + */ + void Init(const TDataset& dataset, index_t start, index_t count, + const EmptyStatistic& left_stat, const EmptyStatistic& right_stat) { } }; #endif diff --git a/fastlib/branches/fastlib-stl/fastlib/tree/tree_test.cc b/fastlib/branches/fastlib-stl/fastlib/tree/tree_test.cc index a174171128..c095e299e1 100644 --- a/fastlib/branches/fastlib-stl/fastlib/tree/tree_test.cc +++ b/fastlib/branches/fastlib-stl/fastlib/tree/tree_test.cc @@ -50,13 +50,13 @@ void TestBallBound() { // Create two balls with a center distance of 1 from each other. // Give the first one a radius of 0.3 and the second a radius of 0.4. - b1.center().Init(3); + b1.center().set_size(3); b1.center()[0] = 1; b1.center()[1] = 2; b1.center()[2] = 3; b1.set_radius(0.3); - b2.center().Init(3); + b2.center().set_size(3); b2.center()[0] = 1; b2.center()[1] = 2; b2.center()[2] = 4; @@ -84,8 +84,7 @@ void TestBallBound() { TEST_ASSERT(!b1.Contains(b2.center())); TEST_ASSERT(!b2.Contains(b1.center())); TEST_ASSERT(b2.Contains(b2.center())); - Vector b2point; // a point that's within the radius bot not the center - b2point.Init(3); + arma::vec b2point(3); // a point that's within the radius bot not the center b2point[0] = 1.1; b2point[1] = 2.1; b2point[2] = 4.1;