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

hi

Browse Source
This commit is contained in:
Garry Boyer
2007-06-03 05:25:15 +00:00
parent f0c70b6108
commit 2e1ed96d8a
8 changed files with 497 additions and 225 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/u/garryb/nbr/affinity.cc
+367 -165
View File
@@ -1,8 +1,17 @@
#include "spbounds.h"
#include "gnp.h"
#include "dfs.h"
#include "nbr_utils.h"
#include "fastlib/fastlib.h"
struct AffinityCommon {
/** The bounding type. Required by NBR. */
typedef SpHrectBound<2> Bound;
/**
* Alpha corresponds to "maximum availability" with the != k condition.
*/
struct Alpha {
public:
double max1;
@@ -16,7 +25,7 @@ struct AffinityCommon {
}
public:
double get(index_t i) {
double get(index_t i) const {
if (unlikely(i == max1_index)) {
return max2;
} else {
@@ -25,24 +34,19 @@ struct AffinityCommon {
}
};
struct AlphaInfo {
struct CombinedInfo {
/** Maximum availability of the point. */
Alpha alpha;
OT_DEF(AlphaInfo) {
OT_MY_OBJECT(alpha);
}
};
struct RhoInfo {
/** Sum of responsibilities. */
double rho;
OT_DEF(RhoInfo) {
OT_DEF(CombinedInfo) {
OT_MY_OBJECT(alpha);
OT_MY_OBJECT(rho);
}
};
typedef SpVectorPoint<AlphaInfo> AlphaPoint;
typedef SpVectorPoint<RhoInfo> RhoPoint;
typedef SpVectorInfoPoint<CombinedInfo> CombinedPoint;
struct Param {
public:
@@ -52,130 +56,144 @@ struct AffinityCommon {
index_t dim;
/** Number of points */
index_t n_points;
/** Self-pereference. */
double pref;
/** The damping factor. */
double lambda;
OT_DEF(Param) {
OT_MY_OBJECT(eps);
OT_MY_OBJECT(dim);
OT_MY_OBJECT(n_points);
OT_MY_OBJECT(pref);
OT_MY_OBJECT(lambda);
}
public:
void Copy(const Param& other) {
pref = other.pref;
eps = other.eps;
dim = other.dim;
n_points = other.n_points;
lambda = other.lambda;
}
void Init(datanode *module) {
dim = -1;
eps = fx_param_double(module, "eps", 1.0e-2);
pref = fx_param_double_req(module, "pref");
lambda = fx_param_double(module, "lambda", 0.8);
}
void AnalyzePoint(const AlphaPoint& q) {
if (dim == -1) {
dim = q.vec().length();
} else {
DEBUG_ASSERT_MSG(dim == q.length(), "Differing dimensionality");
}
void BootstrapMonochromatic(CombinedPoint *point, index_t count) {
dim = point->vec().length();
n_points = count;
// TODO: Realistic values
point->info().rho = 0;
point->info().alpha.max1 = 0;
point->info().alpha.max2 = 0;
point->info().alpha.max1_index = 0;
}
void AnalyzePoint(const RhoPoint& r) {
if (dim == -1) {
dim = q.vec().length();
} else {
DEBUG_ASSERT_MSG(dim == q.length(), "Differing dimensionality");
}
};
struct CombinedStat {
public:
SpRange alpha;
SpRange rho;
OT_DEF(CombinedStat) {
OT_MY_OBJECT(alpha);
OT_MY_OBJECT(rho);
}
public:
double Similarity(const Vector& q, const Vector& r) const {
return 1.0 / sqrt(la::EuclideanDistanceSq(a, b));
void Init(const Param& param) {
Reset(param);
}
double Similarity(
void Reset(const Param& param) {
alpha.InitEmptySet();
rho.InitEmptySet();
}
void Accumulate(const Param& param, const CombinedPoint& point) {
alpha |= SpRange(point.info().alpha.max2, point.info().alpha.max1);
rho |= point.info().rho;
}
void Accumulate(const Param& param,
const CombinedStat& stat, const Bound& bound, index_t n) {
alpha |= stat.alpha;
rho |= stat.rho;
}
void Postprocess(const Param& param, const Bound& bound, index_t n) {}
};
typedef SpNode<Bound, CombinedStat> CombinedNode;
struct Helpers {
static double Similarity(double distsq) {
return -distsq;
}
static double Similarity(const Vector& a, const Vector& b) {
//uint32 anum = (mem::PointerAbsoluteAddress(a.ptr()) * 315187727);
//uint32 bnum = (mem::PointerAbsoluteAddress(b.ptr()) * 210787727);
//uint32 val = ((anum ^ bnum) >> 16) & 0xfff;
//double noise = (1.0e-6 / 4096) * val;
return Similarity(la::DistanceSqEuclidean(a, b));
}
static double Similarity(
const Param& param,
const Vector& q, index_t q_index,
const Vector& r, index_t r_index) const {
const Vector& r, index_t r_index) {
if (unlikely(q_index == r_index)) {
return pref;
return param.pref;
} else {
return Similarity(q, r);
}
return Similarity(q, r);
}
double SimilarityHi(const AlphaNode& a, const AlphaNode& b) const {
double dist = sqrt(a->bound().MinDistanceSqToBound(b->bound()));
double hi = 1.0 / dist;
if (q->begin() < r->end() && r->begin() < q->end() && pref > upper_bound) {
hi = pref;
static double SimilarityHi(
const Param& param,
const CombinedNode& a, const CombinedNode& b) {
double distsq = a.bound().MinDistanceSqToBound(b.bound());
double hi = Similarity(distsq);
if (a.begin() < b.end() && b.begin() < a.end()
&& param.pref > hi) {
hi = param.pref;
}
return hi;
}
static double SimilarityLo(
const Param& param,
const CombinedNode& a, const CombinedNode& b) {
double distsq = a.bound().MaxDistanceSqToBound(b.bound());
double lo = Similarity(distsq);
if (a.begin() < b.end() && b.begin() < a.end()
&& param.pref < lo) {
lo = param.pref;
}
return lo;
}
double ErrorShare(double abs_error_used, const RNode& r_node) {
return (eps - abs_error_used) * r_node.count() / n_points;
static double ErrorShare(const Param& param,
double abs_error_used, const CombinedNode& r_node) {
return (param.eps - abs_error_used) * r_node.count() / param.n_points;
}
};
struct AlphaStat {
public:
SpRange alpha;
OT_DEF(AlphaStat) {
OT_MY_OBJECT(alpha);
}
public:
void Init(const Param& param) {
alpha.InitEmptySet();
}
void Accumulate(const Param& param, const AlphaPoint& point) {
alpha |= SpRange(point.info().max2, point.info().max1);
}
void Accumulate(const Param& param,
const AlphaStat& stat, const Bound& bound, index_t n) {
alpha |= stat.alpha;
}
void Postprocess(const Param& param, const Bound& bound, index_t n) {}
};
struct RhoStat {
public:
SpRange rho;
OT_DEF(RhoStat) {
OT_MY_OBJECT(rho);
}
public:
void Init(const Param& param) {
rhos.InitEmptySet();
}
void Accumulate(const Param& param, const RhoPoint& point) {
rhos |= SpRange(point.info().rho);
}
void Accumulate(const Param& param,
const RhoStat& stat, const Bound& bound, index_t n) {
rhos |= stat.rhos;
}
void Postprocess(const Param& param, const Bound& bound, index_t n) {}
};
typedef SpNode<Bound, AlphaStat> AlphaNode;
typedef SpNode<Bound, RhoStat> RhoNode;
};
class AffinityAlpha {
public:
typedef AffinityCommon::AlphaPoint QPoint;
typedef AffinityCommon::RhoPoint RPoint;
typedef AffinityCommon::CombinedPoint QPoint;
typedef AffinityCommon::CombinedPoint RPoint;
typedef AffinityCommon::Alpha Alpha;
typedef AffinityCommon::Param Param;
typedef AffinityCommon::AlphaStat QStat;
typedef AffinityCommon::RhoStat RStat;
typedef AffinityCommon::AlphaNode QNode;
typedef AffinityCommon::RhoNode RNode;
typedef AffinityCommon::CombinedNode QNode;
typedef AffinityCommon::CombinedNode RNode;
typedef BlankGlobalResult GlobalResult;
struct BlankPostponed QPostponed;
typedef BlankQPostponed QPostponed;
struct Delta {
public:
@@ -200,12 +218,12 @@ class AffinityAlpha {
public:
void Init(const Param& param) {
alpha.max1 = -DBL_MAX;
alpha.max2 = -DBL_MAX;
alpha.max1 = param.pref;
alpha.max2 = param.pref;
alpha.max1_index = -1;
}
void Postprocess(const Param& param,
const QPoint& q, const RNode& r_root) {}
const QPoint& q, index_t q_index, const RNode& r_root) {}
void ApplyPostponed(const Param& param,
const QPostponed& postponed, const QPoint& q) {}
};
@@ -215,9 +233,7 @@ class AffinityAlpha {
SpRange alpha;
OT_DEF(QMassResult) {
OT_MY_OBJECT(alpha_hi);
OT_MY_OBJECT(alpha_lo);
OT_MY_OBJECT(xyz);
OT_MY_OBJECT(alpha);
}
public:
@@ -230,13 +246,13 @@ class AffinityAlpha {
void ApplyDelta(const Param& param, const Delta& delta) {
alpha.MaxWith(delta.alpha);
}
bool ApplyPostponed(const Param& param,
void ApplyPostponed(const Param& param,
const QPostponed& postponed, const QNode& q_node) {}
void StartReaccumulate(const Param& param, const QNode& q_node) {
alpha.InitEmptySet();
}
void Accumulate(const Param& param, const QResult& result) {
alpha |= SpRange(result.max2, result.max1);
alpha |= SpRange(result.alpha.max2, result.alpha.max1);
}
void Accumulate(const Param& param,
const QMassResult& result, index_t n_points) {
@@ -245,6 +261,31 @@ class AffinityAlpha {
void FinishReaccumulate(const Param& param, const QNode& q_node) {}
};
/*
As two-variable functions:
\rho(i, k) = \sum_{j != i, j != k} max(0, S(j,k) - \alpha(j,k))
\alpha(i, k) = min(0, \max_{j != k} S(j,j) + \alpha(j,j) + \rho(i,j))
As one-variable rho:
\rho(k) = \sum_{j != k} max(0, S(j,k) - \alpha(j,k))
\sum max(0, S(j,k) - \alpha(j,k)) - max(0, S(k,k) - alpha(k,k))
\alpha(i, k) = min(0, \max_{j != k} S(j,j) + \alpha(j,j)
+ \rho(i) - max(0, S(i, j) - \alpha(i, j)))
\alpha(i) = min(0, \max^2{j}
S(j,j) + \alpha(j,j) + \rho(j) - max(0, S(i, j) - \alpha(i, j)))
except when i = j in which case we don't need to do the second part
S(j,j) + \alpha(j,j) + \rho(j) - S(i, j) + min(S(i,j), \alpha(i, j))
*/
struct PairVisitor {
public:
Alpha alpha;
@@ -253,32 +294,36 @@ class AffinityAlpha {
void Init(const Param& param) {}
bool StartVisitingQueryPoint(const Param& param,
const Point& q,
const QPoint& q,
const RNode& r_node, const QMassResult& unapplied_mass_results,
QResult* q_result, GlobalResult* global_result) {
alpha = q_result->alpha;
return true;
}
void VisitPair(const Param& param,
const QPoint& q, index_t q_index,
const RPoint& r, index_t r_index) {
double cur_alpha;
if (unlikely(q_index == r_index)) {
cur_alpha = r.info().rho + q.info().alpha.get(r_index);
} else {
double sim = param.Similarity(q.vec(), r.vec());
cur_alpha = min(
const QPoint& q, index_t q_index, const RPoint& r, index_t r_index) {
double candidate_alpha;
if (likely(q_index != r_index)) {
double sim = AffinityCommon::Helpers::Similarity(q.vec(), r.vec());
candidate_alpha = min(
min(sim, q.info().alpha.get(r_index)) + r.info().rho,
sim);
//double availability =
// r.info().rho - math::ClampNonNegative(
// sim - q.info().alpha.get(r_index));
//candidate_alpha = math::ClampNonPositive(availability) + sim;
} else {
candidate_alpha = r.info().rho + q.info().alpha.get(r_index);
}
if (unlikely(cur_alpha > alpha.max2)) {
if (unlikely(cur_alpha > alpha.max1)) {
if (unlikely(candidate_alpha > alpha.max2)) {
if (unlikely(candidate_alpha > alpha.max1)) {
alpha.max2 = alpha.max1;
alpha.max1 = cur_alpha;
alpha.max_index = r_index;
alpha.max1 = candidate_alpha;
alpha.max1_index = r_index;
} else {
alpha.max2 = cur_alpha;
alpha.max2 = candidate_alpha;
}
}
}
@@ -296,14 +341,23 @@ class AffinityAlpha {
const QNode& q_node, const RNode& r_node,
Delta* delta,
GlobalResult* global_result, QPostponed* q_postponed) {
double sim_hi = param.SimilarityHi(q_node, r_node);
double sim_lo = param.SimilarityLo(q_node, r_node);
//WALDO
double sim_lo = AffinityCommon::Helpers::SimilarityLo(
param, q_node, r_node);
delta->alpha.lo = min(
min(q_node.stat().alpha.lo, sim_lo) + r_node.stat().rho.lo,
sim_lo);
delta->alpha.hi =
q_node.stat().alpha.hi + r_node.stat().rho.hi;
if (q_node.begin() < r_node.end() && r_node.begin() < q_node.end()) {
delta->alpha.hi =
q_node.stat().alpha.hi + r_node.stat().rho.hi;
} else {
double sim_hi = AffinityCommon::Helpers::SimilarityHi(
param, q_node, r_node);
delta->alpha.hi = min(
min(q_node.stat().alpha.hi, sim_hi) + r_node.stat().rho.hi,
sim_hi);
}
return true;
}
@@ -325,7 +379,7 @@ class AffinityAlpha {
}
static double Heuristic(const Param& param,
const QNode& q_node, const RNode& r_node, const Delta& delta) {
return -delta.sim_hi;
return -delta.alpha.hi;
}
};
};
@@ -334,16 +388,13 @@ class AffinityRho {
public:
typedef AffinityCommon::Alpha Alpha;
typedef AffinityCommon::RhoPoint QPoint;
typedef AffinityCommon::AlphaPoint RPoint;
typedef AffinityCommon::CombinedPoint QPoint;
typedef AffinityCommon::CombinedPoint RPoint;
typedef AffinityCommon::Param Param;
typedef AffinityCommon::RhoStat QStat;
typedef AffinityCommon::AlphaStat RStat;
typedef AffinityCommon::AlphaNode QNode;
typedef AffinityCommon::RhoNode RNode;
typedef AffinityCommon::CombinedNode QNode;
typedef AffinityCommon::CombinedNode RNode;
typedef BlankGlobalResult GlobalResult;
@@ -359,7 +410,7 @@ class AffinityRho {
public:
void Init(const Param& param) {
Reset();
Reset(param);
}
void Reset(const Param& param) {
@@ -402,7 +453,15 @@ class AffinityRho {
abs_error_used = 0;
}
void Postprocess(const Param& param,
const QPoint& q, const RNode& r_root) {}
const QPoint& q, index_t q_index, const RNode& r_root) {
double responsibility =
param.pref - q.info().alpha.get(q_index);
// Make sure we count ourselves regardless of sign.
if (responsibility < 0) {
rho += responsibility;
}
}
void ApplyPostponed(const Param& param,
const QPostponed& postponed, const QPoint& q) {
rho += postponed.d_rho;
@@ -432,7 +491,7 @@ class AffinityRho {
void ApplyDelta(const Param& param, const Delta& delta) {
rho += delta.d_rho;
}
bool ApplyPostponed(const Param& param,
void ApplyPostponed(const Param& param,
const QPostponed& postponed, const QNode& q_node) {
rho += postponed.d_rho;
abs_error_used += postponed.abs_error_used;
@@ -464,18 +523,19 @@ class AffinityRho {
const RNode& r_node, const QMassResult& unapplied_mass_results,
QResult* q_result, GlobalResult* global_result) {
rho = q_result->rho;
return true;
}
void VisitPair(const Param& param,
const QPoint& q, index_t q_index,
const RPoint& r, index_t r_index) {
double sim = param.Similarity(q.vec(), r.vec())
double responsibility =
AffinityCommon::Helpers::Similarity(
param, q.vec(), q_index, r.vec(), r_index)
- r.info().alpha.get(r_index);
if (sim < 0 && likely(q_index != r_index)) {
sim = 0;
if (responsibility > 0) {
rho += responsibility;
}
rho += sim;
}
void FinishVisitingQueryPoint(const Param& param,
const QPoint& q,
@@ -491,28 +551,42 @@ class AffinityRho {
const QNode& q_node, const RNode& r_node,
Delta* delta,
GlobalResult* global_result, QPostponed* q_postponed) {
double sim_hi = param.SimilarityHi(q_node, r_node);
double sim_lo = param.SimilarityLo(q_node, r_node);
double sim_hi = AffinityCommon::Helpers::SimilarityHi(
param, q_node, r_node);
double sim_lo = AffinityCommon::Helpers::SimilarityLo(
param, q_node, r_node);
delta->rho.lo = (sim_lo - r_node.stat().alpha.hi) * r_node.count();
delta->rho.hi = max(0, sim_hi - r_node.stat().alpha.lo) * r_node.count();
// fprintf(stderr, "(%d,%d) alpha.lo,hi = (%f, %f)\n",
// r_node.begin(), r_node.end(),
// r_node.stat().alpha.hi,
// r_node.stat().alpha.lo
// );
delta->d_rho.lo = max(0.0, sim_lo - r_node.stat().alpha.hi)
* r_node.count();
delta->d_rho.hi = max(0.0, sim_hi - r_node.stat().alpha.lo)
* r_node.count();
return true;
return delta->d_rho.hi != 0;
}
static bool ConsiderPairExtrinsic(const Param& param,
const QNode& q_node, const RNode& r_node, const Delta& delta,
const QMassResult& q_mass_result, const GlobalResult& global_result,
QPostponed* q_postponed) {
double abs_error = delta->d_rho.width() / 2;
/*
double abs_error = delta.d_rho.width() / 2;
double rel_error_hi = abs_error / q_mass_result.rho.lo;
if (rel_error_hi < param.ErrorShare(q_mass_result.abs_error_used, r_node)) {
if (rel_error_hi < AffinityCommon::Helpers::ErrorShare(
param, q_mass_result.abs_error_used, r_node)) {
q_postponed->abs_error_used += abs_error;
q_postponed->d_rho += delta->d_rho.mid();
q_postponed->d_rho += delta.d_rho.mid();
return false;
} else {
return true;
}
*/
return true;
}
static bool ConsiderQueryTermination(const Param& param,
const QNode& q_node,
@@ -523,35 +597,163 @@ class AffinityRho {
static double Heuristic(const Param& param,
const QNode& q_node, const RNode& r_node, const Delta& delta) {
// favor whatever brings our lower bound up the fastest
return -delta.d_rho_lo;
return -delta.d_rho.lo;
}
};
};
void FindExemplars(index_t dimensionality, index_t n_points,
CacheArray<AffinityAlpha::QPoint> *data_points) {
ArrayList<Vector> exemplars;
CacheReadIterator<AffinityAlpha::QPoint> point(data_points, 0);
exemplars.Init();
for (index_t point_i = 0; point_i < n_points; point_i++, point.Next()) {
//ot::Print(point->info());
if (point->info().rho > 0) {
exemplars.AddBack()->Copy(point->vec());
}
}
ot::Print(exemplars);
Matrix m;
m.Init(dimensionality, exemplars.size());
for (index_t i = 0; i < exemplars.size(); i++) {
Vector dest;
m.MakeColumnVector(i, &dest);
dest.CopyValues(exemplars[i]);
}
data::Save("exemplars.txt", m);
}
void AffinityMain(datanode *module, const char *gnp_name) {
typename GNP::Param param;
AffinityAlpha::Param param;
param.Init(fx_submodule(module, gnp_name, gnp_name));
TempCacheArray<typename GNP::QData> q_points;
TempCacheArray<typename GNP::QNode> q_nodes;
TempCacheArray<typename GNP::RData> r_points;
TempCacheArray<typename GNP::RNode> r_nodes;
TempCacheArray<typename GNP::QResult> q_results;
TempCacheArray<AffinityAlpha::QPoint> data_points;
TempCacheArray<AffinityAlpha::QNode> data_nodes;
nbr_utils::LoadKdTree(fx_submodule(module, "q", "q"),
&param, &q_points, &q_nodes);
nbr_utils::LoadKdTree(fx_submodule(module, "r", "r"),
&param, &r_points, &r_nodes);
index_t n_block_points = fx_param_int(
module, "n_block_points", 1024);
index_t n_block_nodes = fx_param_int(
module, "n_block_nodes", 128);
typename GNP::QResult default_result;
default_result.Init(param);
q_results.Init(default_result, q_points.end_index(),
q_points.n_block_elems());
datanode *data_module = fx_submodule(module, "data", "data");
Solver solver;
solver.InitSolve(fx_submodule(module, "solver", "solver"), param, 0,
q_points.cache(), q_nodes.cache(),
r_points.cache(), r_nodes.cache(), q_results.cache());
fx_timer_start(data_module, "read");
Matrix data_matrix;
MUST_PASS(data::Load(fx_param_str_req(data_module, ""), &data_matrix));
index_t n_points = data_matrix.n_cols();
index_t dimensionality = data_matrix.n_rows();
AffinityAlpha::QPoint default_point;
default_point.vec().Init(data_matrix.n_rows());
param.BootstrapMonochromatic(&default_point, data_matrix.n_cols());
data_points.Init(default_point, data_matrix.n_cols(), n_block_points);
for (index_t i = 0; i < data_matrix.n_cols(); i++) {
CacheWrite<AffinityAlpha::QPoint> point(&data_points, i);
point->vec().CopyValues(data_matrix.GetColumnPtr(i));
}
data_matrix.Destruct();
data_matrix.Init(0, 0);
fx_timer_stop(data_module, "read");
AffinityAlpha::QNode data_example_node;
data_example_node.Init(dimensionality, param);
data_nodes.Init(data_example_node, 0, n_block_nodes);
KdTreeMidpointBuilder
<AffinityAlpha::QPoint, AffinityAlpha::QNode, AffinityAlpha::Param>
::Build(data_module, param, &data_points, &data_nodes);
// All the above is not any different for affinity than for anything
// else. Now, time for the iteration.
int n_iter = fx_param_int(module, "n_iter", 1000);
int stable_iterations = 0;
index_t n_changed = n_points;
for (int iter = 0; iter < n_iter && stable_iterations < 50; iter++) {
{
TempCacheArray<AffinityAlpha::QResult> q_results_alpha;
AffinityAlpha::QResult default_result_alpha;
default_result_alpha.Init(param);
q_results_alpha.Init(default_result_alpha, data_points.end_index(),
data_points.n_block_elems());
nbr_utils::ThreadedDualTreeSolver
< AffinityAlpha, DualTreeDepthFirst<AffinityAlpha> >::Solve(
fx_submodule(module, "threads", "iter%d_alpha", iter), param,
&data_points, &data_nodes, &data_points, &data_nodes,
&q_results_alpha);
for (index_t i = 0; i < n_points; i++) {
CacheWrite<AffinityAlpha::QPoint> point(&data_points, i);
CacheRead<AffinityAlpha::QResult> alpha(&q_results_alpha, i);
point->info().alpha = alpha->alpha;
}
nbr_utils::StatFixer<
AffinityAlpha::Param, AffinityAlpha::QPoint, AffinityAlpha::QNode>
::Fix(param, &data_points, &data_nodes);
}
{
TempCacheArray<AffinityRho::QResult> q_results_rho;
AffinityRho::QResult default_result_rho;
default_result_rho.Init(param);
q_results_rho.Init(default_result_rho, data_points.end_index(),
data_points.n_block_elems());
nbr_utils::ThreadedDualTreeSolver
< AffinityRho, DualTreeDepthFirst<AffinityRho> >::Solve(
fx_submodule(module, "threads", "iter%d_rho", iter), param,
&data_points, &data_nodes, &data_points, &data_nodes,
&q_results_rho);
double temperature = 0;//1.0 * n_changed / n_points;
double lambda = param.lambda * (1.0 - temperature) + 0.6 * temperature;
double nonlambda = 1.0 - lambda;
n_changed = 0;
for (index_t i = 0; i < n_points; i++) {
CacheWrite<AffinityAlpha::QPoint> point(&data_points, i);
CacheRead<AffinityRho::QResult> compute_rho(&q_results_rho, i);
double old_rho = point->info().rho;
point->info().rho = old_rho*lambda + compute_rho->rho*nonlambda;
if ((old_rho > 0) != (point->info().rho > 0)) {
n_changed++;
}
}
if (n_changed == 0) {
stable_iterations++;
} else {
stable_iterations = 0;
}
fprintf(stderr, "------------- Changed = %"LI"d\n", n_changed);
nbr_utils::StatFixer<
AffinityAlpha::Param, AffinityAlpha::QPoint, AffinityAlpha::QNode>
::Fix(param, &data_points, &data_nodes);
}
}
FindExemplars(dimensionality, n_points, &data_points);
}
int main(int argc, char *argv[]) {
fx_init(argc, argv);
AffinityMain(fx_root, "affinity");
fx_done();
}
fastlib/u/garryb/nbr/allnn.cc
+1 -1
View File
@@ -84,7 +84,7 @@ class Allnn {
}
void Postprocess(const Param& param,
const QPoint& q_point,
const QPoint& q_point, index_t q_index,
const RNode& r_root) {}
void ApplyPostponed(const Param& param,
const QPostponed& postponed,
fastlib/u/garryb/nbr/build.py
+4
View File
@@ -32,6 +32,10 @@ binrule(name = "allnn",
sources = ["allnn.cc"],
deplibs = [":nbr"])
binrule(name = "affinity",
sources = ["affinity.cc"],
deplibs = [":nbr"])
binrule(name = "cache_test",
fastlib/u/garryb/nbr/cache.h
+2 -2
View File
@@ -352,13 +352,13 @@ class CacheArray {
}
void ReleaseBlock_(BlockDevice::blockid_t blockid) {
DEBUG_ONLY(--metadatas_[blockid].lock_count);
DEBUG_ONLY(--metadatas_[blockid - begin_block_fake_].lock_count);
}
void ReleaseElement_(index_t element_id) {
DEBUG_ONLY(BoundsCheck_(element_id));
DEBUG_ONLY(
ReleaseBlock_((element_id >> n_block_elems_log_) - begin_block_fake_));
ReleaseBlock_((element_id >> n_block_elems_log_)));
}
};
fastlib/u/garryb/nbr/dfs.h
+5 -6
View File
@@ -185,7 +185,7 @@ void DualTreeDepthFirst<GNP>::PushDown_(
const typename GNP::QPoint *q_point = q_points_.StartRead(q_i);
q_result->ApplyPostponed(param_, q_node_mut->postponed, *q_point);
q_result->Postprocess(param_, *q_point, *r_root_);
q_result->Postprocess(param_, *q_point, q_i, *r_root_);
q_results_.StopWrite(q_i);
q_points_.StopRead(q_i);
}
@@ -356,13 +356,12 @@ void DualTreeDepthFirst<GNP>::BaseCase_(
q_result->ApplyPostponed(param_, q_node_mut->postponed, *q_point);
if (unlikely(
visitor.StartVisitingQueryPoint(param_, *q_point, *r_node,
exclusive_unvisited, q_result, &global_result_))) {
if (visitor.StartVisitingQueryPoint(param_, *q_point, *r_node,
exclusive_unvisited, q_result, &global_result_)) {
CacheReadIterator<typename GNP::RPoint> r_iter(&r_points_, r_node->begin());
index_t r_end = r_node->end();
for (index_t r_i = r_node->begin(); r_i < r_end; ++r_i, r_iter.Next()) {
for (index_t r_i = r_node->begin(); r_i < r_node->end(); ++r_i,
r_iter.Next()) {
const typename GNP::RPoint *r_point = r_iter;
visitor.VisitPair(param_, *q_point, q_i, *r_point, r_i);
fastlib/u/garryb/nbr/gnp.h
+5 -5
View File
@@ -47,15 +47,15 @@ struct BlankQResult {
public:
OT_DEF(BlankQResult) {}
public:
template<typename Param, typename Point, typename QPointInfo, typename RNode>
template<typename Param, typename Point, typename RNode>
void Init(const Param& param,
const Point& q_point, const QPointInfo& q_info,
const Point& q_point,
const RNode& r_root) {}
template<typename Param, typename Point, typename QPointInfo, typename RNode>
template<typename Param, typename Point, typename RNode>
void Postprocess(const Param& param,
const Point& q_point, const QPointInfo& q_info,
const Point& q_point,
const RNode& r_root) {}
template<typename Param, typename QPostponed, typename QPointInfo, typename Point>
template<typename Param, typename QPostponed, typename Point>
void ApplyPostponed(const Param& param,
const QPostponed& postponed,
const Point& q_point) {}
fastlib/u/garryb/nbr/nbr_utils.h
+64 -2
View File
@@ -10,6 +10,68 @@
namespace nbr_utils {
template<typename Param, typename Point, typename Node>
class StatFixer {
public:
static void Fix(const Param &param,
CacheArray<Point> *points, CacheArray<Node> *nodes) {
StatFixer fixer;
fixer.InitFix(&param, points, nodes);
}
private:
const Param *param_;
CacheArray<Point> points_;
CacheArray<Node> nodes_;
public:
void InitFix(const Param *param,
CacheArray<Point> *points, CacheArray<Node> *nodes);
private:
void FixRecursively_(index_t node_index);
};
template<typename Param, typename Point, typename Node>
void StatFixer<Param, Point, Node>::InitFix(
const Param *param, CacheArray<Point> *points, CacheArray<Node> *nodes) {
param_ = param;
points_.Init(points, BlockDevice::READ);
nodes_.Init(nodes, BlockDevice::MODIFY);
FixRecursively_(0);
nodes_.Flush();
points_.Flush();
}
template<typename Param, typename Point, typename Node>
void StatFixer<Param, Point, Node>::FixRecursively_(index_t node_index) {
CacheWrite<Node> node(&nodes_, node_index);
node->stat().Reset(*param_);
if (!node->is_leaf()) {
for (index_t k = 0; k < 2; k++) {
index_t child_index = node->child(k);
FixRecursively_(child_index);
CacheRead<Node> child(&nodes_, child_index);
node->stat().Accumulate(*param_, child->stat(),
child->bound(), child->count());
}
node->stat().Postprocess(*param_, node->bound(),
node->count());
} else {
CacheReadIterator<Point> point(&points_, node->begin());
for (index_t i = 0; i < node->count(); i++, point.Next()) {
node->stat().Accumulate(*param_, *point);
}
}
node->stat().Postprocess(*param_, node->bound(), node->count());
}
template<typename GNP, typename Solver>
class ThreadedDualTreeSolver {
private:
@@ -42,7 +104,7 @@ class ThreadedDualTreeSolver {
name.InitSprintf("grain_%d", work[i]);
base_->mutex_.Lock();
struct datanode *submodule = fx_submodule(base_->module_,
name.c_str(), name.c_str());
name.c_str(), "solver");
base_->mutex_.Unlock();
solver.InitSolve(submodule, *base_->param_, q_root_index,
@@ -213,7 +275,7 @@ struct MpiDualTreeConfig {
void Copy(const MpiDualTreeConfig& other) {
*this = other;
KK }
}
OT_DEF(MpiDualTreeConfig) {
OT_MY_OBJECT(n_threads);
fastlib/u/garryb/nbr/tkde.cc
+49 -44
View File
@@ -41,18 +41,14 @@ class Tkde {
public:
/** The bounding type. Required by NBR. */
typedef SpHrectBound<2> Bound;
/** The type of point in use. Required by NBR. */
typedef Vector Point;
typedef SpVectorPoint QPoint;
typedef SpVectorPoint RPoint;
/** The type of kernel in use. NOT required by NBR. */
typedef EpanKernel Kernel;
/** Per-query statistic. Required by NBR. */
typedef BlankStat QStat;
typedef BlankGlobalResult GlobalResult;
typedef BlankPointInfo QPointInfo;
typedef BlankPointInfo RPointInfo;
/**
* All parameters required by the execution of the algorithm.
@@ -62,21 +58,25 @@ class Tkde {
struct Param {
public:
/**
* The threshold in use.
* This is a range to allow for epsilon checking.
* The threshold in use, with upper and lower bounds to prevent
* roundoff error.
*
* This is also normalized for dimensionality and the number of reference
* points.
*/
SpRange thresh;
/** The kernel in use. */
Kernel kernel;
/** The dimensionality of the data sets. */
index_t dim;
/** The original threshold */
double threshold_orig;
/** The original threshold, before normalization. */
double nominal_threshold;
OT_DEF(Param) {
OT_MY_OBJECT(kernel);
OT_MY_OBJECT(thresh);
OT_MY_OBJECT(kernel);
OT_MY_OBJECT(dim);
OT_MY_OBJECT(nominal_threshold);
}
public:
@@ -90,20 +90,16 @@ class Tkde {
* Initialize parameters from a data node (Req NBR).
*/
void Init(datanode *module) {
dim = -1;
kernel.Init(fx_param_double_req(module, "h"));
threshold_orig = fx_param_double_req(module, "threshold");
nominal_threshold = fx_param_double_req(module, "threshold");
}
void AnalyzePoint(const Point& q_point) {
if (dim == -1) {
dim = q_point.length();
double t = threshold_orig * kernel.CalcNormConstant(dim);
thresh.lo = t * (1.0 - 1.0e-4);
thresh.hi = t * (1.0 + 1.0e-4);
} else {
DEBUG_ASSERT_MSG(dim == q_point.length(), "Differing dimensionality");
}
void BootstrapMonochromatic(QPoint* point, index_t count) {
dim = point->vec().length();
double normalized_threshold =
nominal_threshold * kernel.CalcNormConstant(dim) * count;
thresh.lo = normalized_threshold * (1.0 - 1.0e5);
thresh.hi = normalized_threshold * (1.0 + 1.0e5);
}
public:
@@ -114,11 +110,11 @@ class Tkde {
* the actual query point (not NBR).
*/
double ComputeKernelSum(
const Vector& q_point,
const Vector& q,
index_t r_count, const Vector& r_mass, double r_sumsq) const {
double quadratic_term =
+ r_count * la::Dot(q_point, q_point)
- 2.0 * la::Dot(q_point, r_mass)
+ r_count * la::Dot(q, q)
- 2.0 * la::Dot(q, r_mass)
+ r_sumsq;
return r_count - quadratic_term * kernel.inv_bandwidth_sq();
}
@@ -245,8 +241,8 @@ class Tkde {
/**
* Accumulate data from a single point (Req NBR).
*/
void Accumulate(const Param& param, const Vector& point) {
moment_info.Add(1, point, la::Dot(point, point));
void Accumulate(const Param& param, const QPoint& point) {
moment_info.Add(1, point.vec(), la::Dot(point.vec(), point.vec()));
}
/**
@@ -264,8 +260,18 @@ class Tkde {
void Postprocess(const Param& param, const Bound& bound, index_t n) {
}
};
/**
* Query-tree statistic.
*
* Note that this statistic is not actually needed and a blank statistic
* is fine, but QStat must equal RStat in order for us to allow
* monochromatic execution.
*
* This limitation may be removed in a further version of NBR.
*/
typedef RStat QStat;
/**
* Query node.
*/
@@ -350,7 +356,7 @@ class Tkde {
}
void Postprocess(const Param& param,
const Vector& q_point, const QPointInfo& q_info,
const QPoint& q,
const RNode& r_root) {
if (density > param.thresh.hi) {
label &= LAB_HI;
@@ -362,12 +368,13 @@ class Tkde {
void ApplyPostponed(const Param& param,
const QPostponed& postponed,
const Vector& q_point) {
const QPoint& q) {
label &= postponed.label; /* bitwise OR */
DEBUG_ASSERT(label != LAB_NEITHER);
if (!postponed.moment_info.is_empty()) {
density += postponed.moment_info.ComputeKernelSum(param, q_point);
density += postponed.moment_info.ComputeKernelSum(
param, q.vec());
}
}
};
@@ -463,8 +470,7 @@ class Tkde {
// - this function must assume that global_result is incomplete (which is
// reasonable in allnn)
bool StartVisitingQueryPoint(const Param& param,
const Vector& q_point,
const QPointInfo& q_info,
const QPoint& q,
const RNode& r_node,
const QMassResult& unapplied_mass_results,
QResult* q_result,
@@ -473,17 +479,17 @@ class Tkde {
return false;
}
double distance_sq_lo = r_node.bound().MinDistanceSqToPoint(q_point);
double distance_sq_lo = r_node.bound().MinDistanceSqToPoint(q.vec());
if (unlikely(distance_sq_lo > param.kernel.bandwidth_sq())) {
return false;
}
double distance_sq_hi = r_node.bound().MaxDistanceSqToPoint(q_point);
double distance_sq_hi = r_node.bound().MaxDistanceSqToPoint(q.vec());
if (unlikely(distance_sq_hi < param.kernel.bandwidth_sq())) {
q_result->density += r_node.stat().moment_info.ComputeKernelSum(
param, q_point);
param, q.vec());
return false;
}
@@ -493,15 +499,14 @@ class Tkde {
}
void VisitPair(const Param& param,
const Vector& q_point, const QPointInfo& q_info, index_t q_index,
const Vector& r_point, const RPointInfo& r_info, index_t r_index) {
double distance = la::DistanceSqEuclidean(q_point, r_point);
const QPoint& q, index_t q_index,
const RPoint& r, index_t r_index) {
double distance = la::DistanceSqEuclidean(q.vec(), r.vec());
density += param.kernel.EvalUnnormOnSq(distance);
}
void FinishVisitingQueryPoint(const Param& param,
const Vector& q_point,
const QPointInfo& q_info,
const QPoint& q,
const RNode& r_node,
const QMassResult& unapplied_mass_results,
QResult* q_result,
@@ -627,7 +632,7 @@ int main(int argc, char *argv[]) {
fx_root, "tkde");
MPI_Finalize();
#else
nbr_utils::ThreadedDualTreeMain<Tkde, DualTreeDepthFirst<Tkde> >(
nbr_utils::MonochromaticDualTreeMain<Tkde, DualTreeDepthFirst<Tkde> >(
fx_root, "tkde");
#endif