apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge fastlib-armadillo/fastlib/tree with fastlib-stl/fastlib/tree. This should

Browse Source 
have been done earlier to prevent duplicated work (some functions had been
switched to arma::vec and std::vector in different functions; I stuck with
arma::vec because fixed-size vectors made more sense in those particular
cases)
This commit is contained in:
Ryan Curtin
2010-06-18 17:59:54 +00:00
parent 0c216c8c7b
commit 7e384edcaa
13 changed files with 1877 additions and 1114 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/branches/fastlib-stl/fastlib/tree/CMakeLists.txt
+6
View File
@@ -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
)
fastlib/branches/fastlib-stl/fastlib/tree/bounds.h
+4 -874
View File
@@ -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<int t_pow = 2>
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<DRange>(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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(v_lo); // v_lo is non-negative
sum_hi += math::Pow<t_pow, 1>(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<t_pow, 1>((v1 + fabs(v1)) + (v2 + fabs(v2)));
sum_hi += math::Pow<t_pow, 1>(-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:
* <code>
* other.CalcMidpoint(&other_midpoint)
* return MinDistanceSqToPoint(other_midpoint)
* </code>
*/
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<t_pow, 1>(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<t_pow, 1>(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<t_pow, 1>(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<int t_pow>
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<t_pow>(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<int t_result_pow>
static double PowDistance(const Vector& a, const Vector& b) {
return math::Pow<t_result_pow, t_pow>(
la::RawLMetric<t_pow>(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 TMetric = LMetric<2>, 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:
* <code>
* other.CalcMidpoint(&other_midpoint)
* return MinDistanceSqToPoint(other_midpoint)
* </code>
*/
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
fastlib/branches/fastlib-stl/fastlib/tree/dballbound.h
+172
View File
@@ -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 <armadillo>
/**
* 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 TMetric = LMetric<2>, 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:
* <code>
* other.CalcMidpoint(&other_midpoint)
* return MinDistanceSqToPoint(other_midpoint)
* </code>
*/
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
fastlib/branches/fastlib-stl/fastlib/tree/dballbound_impl.h
+209
View File
@@ -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 <armadillo>
#include "../base/arma_compat.h"
/**
* Determines if a point is within the bound.
*/
template<typename TMetric, typename TPoint>
bool DBallBound<TMetric, TPoint>::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<typename TMetric, typename TPoint>
void DBallBound<TMetric, TPoint>::CalculateMidpoint(Point *centroid) const {
ot::InitCopy(centroid, center_);
}
/**
* Calculates minimum bound-to-point squared distance.
*/
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MinDistance(const Point& point) const {
return math::ClampNonNegative(MidDistance(point) - radius_);
}
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MinDistanceSq(const Point& point) const {
return math::Pow<2, 1>(MinDistance(point));
}
/**
* Calculates minimum bound-to-bound squared distance.
*/
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MinDistance(const DBallBound& other) const {
double delta = MidDistance(other.center_) - radius_ - other.radius_;
return math::ClampNonNegative(delta);
}
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MinDistanceSq(const DBallBound& other) const {
return math::Pow<2, 1>(MinDistance(other));
}
/**
* Computes maximum distance.
*/
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MaxDistance(const Point& point) const {
return MidDistance(point) + radius_;
}
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MaxDistanceSq(const Point& point) const {
return math::Pow<2, 1>(MaxDistance(point));
}
/**
* Computes maximum distance.
*/
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MaxDistance(const DBallBound& other) const {
return MidDistance(other.center_) + radius_ + other.radius_;
}
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::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<typename TMetric, typename TPoint>
DRange DBallBound<TMetric, TPoint>::RangeDistance(const DBallBound& other) const {
double delta = MidDistance(other.center_);
double sumradius = radius_ + other.radius_;
return DRange(
math::ClampNonNegative(delta - sumradius),
delta + sumradius);
}
template<typename TMetric, typename TPoint>
DRange DBallBound<TMetric, TPoint>::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:
* <code>
* other.CalcMidpoint(&other_midpoint)
* return MinDistanceSqToPoint(other_midpoint)
* </code>
*/
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MinToMid(const DBallBound& other) const {
double delta = MidDistance(other.center_) - radius_;
return math::ClampNonNegative(delta);
}
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::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<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MinimaxDistance(const DBallBound& other) const {
double delta = MidDistance(other.center_) + other.radius_ - radius_;
return math::ClampNonNegative(delta);
}
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MinimaxDistanceSq(const DBallBound& other) const {
return math::Pow<2, 1>(MinimaxDistance(other));
}
/**
* Calculates midpoint-to-midpoint bounding box distance.
*/
template< >
inline double DBallBound<LMetric<2>, GenVector<double> >::MidDistance(const GenVector<double>& point) const {
arma::vec tmp1;
arma::vec tmp2;
arma_compat::vectorToVec(center_, tmp1);
arma_compat::vectorToVec(point, tmp2);
return Metric::Distance(tmp1, tmp2);
}
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MidDistance(const DBallBound& other) const {
return MidDistance(other.center_);
}
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MidDistanceSq(const DBallBound& other) const {
return math::Pow<2, 1>(MidDistance(other));
}
template<typename TMetric, typename TPoint>
double DBallBound<TMetric, TPoint>::MidDistance(const Point& point) const {
return Metric::Distance(center_, point);
}
#endif
fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound.h
+234
View File
@@ -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 <armadillo>
/**
* Hyper-rectangle bound for an L-metric.
*
* Template parameter t_pow is the metric to use; use 2 for Euclidean (L2).
*/
template<int t_pow = 2>
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:
* <code>
* other.CalcMidpoint(&other_midpoint)
* return MinDistanceSqToPoint(other_midpoint)
* </code>
*/
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
fastlib/branches/fastlib-stl/fastlib/tree/dhrectbound_impl.h
+687
View File
@@ -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 <math.h>
using arma::vec;
/**
* Empty constructor
*/
template<int t_pow>
DHrectBound<t_pow>::DHrectBound() {
bounds_ = NULL;
dim_ = 0;
}
/**
* Initializes to specified dimensionality with each dimension the empty
* set.
*/
template<int t_pow>
DHrectBound<t_pow>::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<int t_pow>
DHrectBound<t_pow>::~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<int t_pow>
void DHrectBound<t_pow>::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<int t_pow>
void DHrectBound<t_pow>::Reset() {
for (index_t i = 0; i < dim_; i++) {
bounds_[i].InitEmptySet();
}
}
/**
* Sets the dimensionality.
*/
template<int t_pow>
void DHrectBound<t_pow>::SetSize(index_t dim) {
if(bounds_)
delete[] bounds_;
bounds_ = new DRange[dim];
dim_ = dim;
Reset();
}
/**
* Determines if a point is within this bound.
*/
template<int t_pow>
bool DHrectBound<t_pow>::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<int t_pow>
const DRange DHrectBound<t_pow>::operator[](index_t i) const {
return bounds_[i];
}
/**
* Calculates the maximum distance within the rectangle
*/
template<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
void DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<t_pow, 1>(v); // v is non-negative
}
return math::Pow<2, t_pow>(sum);
}
*/
/**
* Computes maximum distance.
*/
template<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
DRange DHrectBound<t_pow>::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<int t_pow>
DRange DHrectBound<t_pow>::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:
* <code>
* other.CalcMidpoint(&other_midpoint)
* return MinDistanceSqToPoint(other_midpoint)
* </code>
*/
template<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
double DHrectBound<t_pow>::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<int t_pow>
DHrectBound<t_pow>& DHrectBound<t_pow>::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<int t_pow>
DHrectBound<t_pow>& DHrectBound<t_pow>::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<int t_pow>
DHrectBound<t_pow>& DHrectBound<t_pow>::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<int t_pow>
DHrectBound<t_pow>& DHrectBound<t_pow>::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
fastlib/branches/fastlib-stl/fastlib/tree/kdtree.h
+138 -103
View File
@@ -47,15 +47,14 @@
#include "spacetree.h"
#include "bounds.h"
//#include "spacetree.h"
//#include "bounds.h"
#include <vector>
#include "../fx/fx.h"
#include "kdtree_impl.h"
//#include "kdtree_impl.h"
#include <armadillo>
/**
* 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<typename TKdTree, typename T>
TKdTree *MakeKdTreeMidpointSelective(GenMatrix<T>& matrix,
const Vector& split_dimensions,
index_t leaf_size,
std::vector<index_t> *old_from_new = NULL,
std::vector<index_t> *new_from_old = NULL) {
TKdTree *MakeKdTreeMidpointSelective(arma::Mat<T>& matrix,
const arma::uvec& split_dimensions,
index_t leaf_size,
arma::Col<index_t>& old_from_new,
arma::Col<index_t>& 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<typename TKdTree, typename T>
TKdTree *MakeKdTreeMidpointSelective(arma::Mat<T>& matrix,
const arma::uvec& split_dimensions,
index_t leaf_size,
arma::Col<index_t>& 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<typename TKdTree, typename T>
TKdTree *MakeKdTreeMidpointSelective(arma::Mat<T>& 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<typename TKdTree, typename T>
TKdTree *MakeKdTreeMidpointSelective(GenMatrix<T>& matrix,
const Vector& split_dimensions,
index_t leaf_size,
GenVector<index_t> *old_from_new = NULL,
GenVector<index_t> *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<typename TKdTree, typename T>
TKdTree *MakeKdTreeMidpoint(GenMatrix<T>& matrix, index_t leaf_size,
std::vector<index_t> *old_from_new = NULL,
std::vector<index_t> *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<T>& matrix,
index_t leaf_size,
arma::Col<index_t>& old_from_new,
arma::Col<index_t>& 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<TKdTree>(matrix, split_dimensions,
leaf_size, old_from_new, new_from_old);
result = MakeKdTreeMidpointSelective<TKdTree>(matrix,
split_dimensions, leaf_size, old_from_new, new_from_old);
return result;
}
template<typename TKdTree, typename T>
TKdTree *MakeKdTreeMidpoint(GenMatrix<T>& matrix,
index_t leaf_size,
GenVector<index_t> *old_from_new = NULL,
GenVector<index_t> *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<T>& matrix,
index_t leaf_size,
arma::Col<index_t>& 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<TKdTree>(matrix,
split_dimensions,
leaf_size, old_from_new, new_from_old);
split_dimensions, leaf_size, old_from_new);
return result;
}
template<typename TKdTree, typename T>
TKdTree *MakeKdTreeMidpoint(arma::Mat<T>& 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<TKdTree>(matrix,
split_dimensions, leaf_size);
return result;
}
@@ -208,7 +242,7 @@ namespace tree {
* @code
* MyKdTree *q_tree;
* Matrix q_matrix;
* std::vector<index_t> q_permutation;
* arma::Col<index_t> 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<typename TKdTree, typename T>
success_t LoadKdTree(datanode *module,
GenMatrix<T> *matrix, TKdTree **tree_pp,
std::vector<index_t> *old_from_new) {
success_t LoadKdTree(datanode* module,
arma::Mat<T>& matrix,
TKdTree** tree_pp,
arma::Col<index_t>& 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<TKdTree>(
*matrix, leaflen, old_from_new);
matrix, leaflen, old_from_new);
fx_timer_stop(module, "make_tree");
}
fx_timer_stop(module, "load");
fastlib/branches/fastlib-stl/fastlib/tree/kdtree_impl.cc
+21
View File
@@ -0,0 +1,21 @@
/**
* @file kdtree_impl.cc
*
* Specialized functions.
*/
#include "../fastlib.h"
#include "../la/matrix.h"
#include "kdtree.h"
#include <armadillo>
/***
* 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]];
}
fastlib/branches/fastlib-stl/fastlib/tree/kdtree_impl.h
+281 -115
View File
@@ -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<typename T>
void MakeBoundVector(const GenVector<T>& point,
const Vector& bound_dimensions, GenVector<T>* 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<typename TBound, typename T>
void SelectFindBoundFromMatrix(const GenMatrix<T>& matrix,
const Vector& split_dimensions, index_t first, index_t count,
TBound *bounds){
void SelectFindBoundFromMatrix(const arma::Mat<T>& 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<T> 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<T> 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<typename TBound, typename T>
void FindBoundFromMatrix(const GenMatrix<T>& 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<T>& 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<typename TBound>
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<typename TBound>
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<index_t>& 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<typename TBound>
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<typename TBound, typename T>
index_t SelectMatrixPartition(GenMatrix<T>& 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<T>& 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<T> 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<T> 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<T> 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<T> 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<T> left_vector;
GenVector<T> 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<T> 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<T> 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<typename TBound, typename T>
index_t SelectMatrixPartition(arma::Mat<T>& 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<index_t>& 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<typename TKdTree, typename T>
void SplitKdTreeMidpoint(GenMatrix<T>& matrix,
TKdTree *node, index_t leaf_size, index_t *old_from_new){
void SplitKdTreeMidpoint(arma::Mat<T>& matrix,
TKdTree *node,
index_t leaf_size,
arma::Col<T>& 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<typename TKdTree, typename T>
void SplitKdTreeMidpoint(arma::Mat<T>& 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<typename TKdTree, typename T>
void SelectSplitKdTreeMidpoint(GenMatrix<T>& matrix,
const Vector& split_dimensions, TKdTree *node, index_t leaf_size,
index_t *old_from_new) {
void SelectSplitKdTreeMidpoint(arma::Mat<T>& matrix,
const arma::uvec& split_dimensions,
TKdTree *node,
index_t leaf_size,
arma::Col<index_t>& 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<typename TKdTree, typename T>
void SelectSplitKdTreeMidpoint(arma::Mat<T>& 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
fastlib/branches/fastlib-stl/fastlib/tree/lmetric.h
+83
View File
@@ -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 <math.h>
#include <armadillo>
/**
* 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<int t_pow>
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<int t_result_pow>
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
fastlib/branches/fastlib-stl/fastlib/tree/spacetree.h
+25 -2
View File
@@ -42,7 +42,9 @@
#include "../base/base.h"
#include "statistic.h"
//#include "statistic.h"
#include <armadillo>
#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_;
}
fastlib/branches/fastlib-stl/fastlib/tree/statistic.h
+14 -16
View File
@@ -50,24 +50,22 @@
*/
template<class TDataset>
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
fastlib/branches/fastlib-stl/fastlib/tree/tree_test.cc
+3 -4
View File
@@ -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;