diff --git a/fastlib/trunk/contrib/tqlong/matching/CMakeLists.txt b/fastlib/trunk/contrib/tqlong/matching/CMakeLists.txt index 9997b11ce4..7f4ad0efe5 100644 --- a/fastlib/trunk/contrib/tqlong/matching/CMakeLists.txt +++ b/fastlib/trunk/contrib/tqlong/matching/CMakeLists.txt @@ -4,8 +4,10 @@ set(SOURCES # auction_max_weight_matching # naive_distance_matrix # kdtree_distance_matrix - matching.cc - kdnode.cc + matching + kdnode + single_tree + auction_matching ) # add directory name to sources diff --git a/fastlib/trunk/contrib/tqlong/matching/auction_matching.h b/fastlib/trunk/contrib/tqlong/matching/auction_matching.h new file mode 100644 index 0000000000..75c145e93d --- /dev/null +++ b/fastlib/trunk/contrib/tqlong/matching/auction_matching.h @@ -0,0 +1,120 @@ +#ifndef AUCTION_MATCHING_H +#define AUCTION_MATCHING_H + +#include +#include "matching.h" + +MATCHING_NAMESPACE_BEGIN; + +template + class AuctionMatching +{ +public: + typedef W weight_type; +protected: + weight_type& weight_; + double epsilon_; + std::vector left_, right_, bidders_; + std::vector bids_; + + void clearMatches(); + void forwardAuction(double &pruned, double &total); + void placeBid(int l, int r, double price); + void setMatch(int l, int r, double price); +public: + AuctionMatching(weight_type& weight); + void doMatch(); + int n_left() const { return weight_.n_rows(); } + int n_right() const { return weight_.n_cols(); } + int left(int index) const { return left_.at(index); } + int right(int index) const { return right_.at(index); } + double getP(int l, int r) { return weight_.get(l, r)-weight_.price(r); } + double price(int r) { return weight_.price(r); } +}; + +template + AuctionMatching::AuctionMatching(weight_type &weight) + : weight_(weight), left_(n_left()), right_(n_right()), + bidders_(n_right()), bids_(n_right()) +{ +} + +template + void AuctionMatching::doMatch() +{ + // epsilon scaling + double total_pruned = 0, total_cals = 0; + for (epsilon_ = 1.0/n_left(); epsilon_ >= 1.0/n_left(); epsilon_ /= 2) + { + double pruned, total; + clearMatches(); + forwardAuction(pruned, total); + total_pruned += pruned; + total_cals += total; + std::cout << "epsilon = " << epsilon_ << " cals = " << total_pruned << "/" << total_cals << "\n"; + } +} + +template + void AuctionMatching::clearMatches() +{ + std::fill(left_.begin(), left_.end(), -1); + std::fill(right_.begin(), right_.end(), -1); +} + +template + void AuctionMatching::forwardAuction(double &pruned, double &total) +{ + pruned = total = 0; + while (1) + { + bool doneMatching = true; + std::fill(bids_.begin(), bids_.end(), -std::numeric_limits::infinity()); + std::fill(bidders_.begin(), bidders_.end(), -1); + weight_.refresh(); + for (int i = 0; i < n_left(); i++) if (left(i) == -1) + { + doneMatching = false; + std::vector bests(2, -1); + pruned += weight_.kBest(i, bests); + total += n_right(); + if (bests[0] == -1 || bests[1] == -1) + { + printf("error kBest\n"); + return; + } + double v = getP(i, bests[0]), w = getP(i, bests[1]); + placeBid(i, bests[0], price(bests[0])+v-w+epsilon_); + } + if (doneMatching) break; + for (int j = 0; j < n_right(); j++) if (bidders_[j] != -1) + { + setMatch(bidders_[j], j, bids_[j]); + } + } +} + +template + void AuctionMatching::placeBid(int l, int r, double price) +{ + if (bids_[r] < price) + { + bids_[r] = price; + bidders_[r] = l; + } +} + +template + void AuctionMatching::setMatch(int l, int r, double price) +{ + weight_.setPrice(r, price); + int ol = right(r); + left_[l] = r; + right_[r] = l; + if (ol != -1) + left_[ol] = -1; +} + +MATCHING_NAMESPACE_END; + +#endif // AUCTION_MATCHING_H diff --git a/fastlib/trunk/contrib/tqlong/matching/kdnode.cc b/fastlib/trunk/contrib/tqlong/matching/kdnode.cc index 5c63bffd59..903a4f807d 100644 --- a/fastlib/trunk/contrib/tqlong/matching/kdnode.cc +++ b/fastlib/trunk/contrib/tqlong/matching/kdnode.cc @@ -41,6 +41,12 @@ int KDNode::oldIndex(int index) const return oldIndex_.at(index+dfsIndex_); } +int KDNode::index(int index) const +{ + if (index < 0 || index >= n_points_) return -1; // error + return index+dfsIndex_; +} + void KDNode::getPoint(int index, Vector &point) const { points_.MakeColumnVector(oldIndex(index), &point); diff --git a/fastlib/trunk/contrib/tqlong/matching/kdnode.h b/fastlib/trunk/contrib/tqlong/matching/kdnode.h index 51226ccf47..4572e13659 100644 --- a/fastlib/trunk/contrib/tqlong/matching/kdnode.h +++ b/fastlib/trunk/contrib/tqlong/matching/kdnode.h @@ -41,6 +41,7 @@ public: int n_points() const; int n_dim() const; + int index(int idx) const; int oldIndex(int index) const; void getPoint(int index, Vector& point) const; double get(int dim, int index) const; @@ -92,63 +93,8 @@ public: std::string toString(int depth = 0) const; }; -template - KDNodeStats::KDNodeStats(const Matrix &points) - : KDNode(points), pointStats_(*(new all_point_stats_type(n_points_))), changed_(true) -{ -} - -template - KDNodeStats::KDNodeStats(KDNodeStats *parent) - : KDNode(parent), pointStats_(parent->pointStats_), changed_(true) -{ -} - -template - KDNodeStats::~KDNodeStats() -{ - if (isRoot()) delete &pointStats_; -} - -template - KDNode* KDNodeStats::newNode(KDNode *parent) -{ - return new KDNodeStats((KDNodeStats*) parent); -} - -template - void KDNodeStats::setPointStats(int index, const point_stats_type &stats) -{ - pointStats_[oldIndex(index)] = stats; - ((KDNodeStats*)leaf(index))->setChanged(true); -} - -template - void KDNodeStats::setChanged(bool changed) -{ - if (changed_ == changed) return; - changed_ = changed; - if (changed && !isRoot()) - ((KDNodeStats*)parent())->setChanged(changed); -} - -template - void KDNodeStats::visit(bool init) -{ - if (!isChanged()) return; // the node is unchanged, not neccessary to proceed - if (isLeaf()) - setLeafStats(init); - else - { - for (unsigned int i = 0; i < children_.size(); i++) - { - ((KDNodeStats*) children_[i])->visit(init); - } - setNonLeafStats(init); - } - setChanged(false); // the node statistics is refreshed -} - MATCHING_NAMESPACE_END; +#include "kdnode_impl.h" + #endif // KDNODE_H diff --git a/fastlib/trunk/contrib/tqlong/matching/kdnode_impl.h b/fastlib/trunk/contrib/tqlong/matching/kdnode_impl.h new file mode 100644 index 0000000000..9186971602 --- /dev/null +++ b/fastlib/trunk/contrib/tqlong/matching/kdnode_impl.h @@ -0,0 +1,67 @@ +#ifndef KDNODE_IMPL_H +#define KDNODE_IMPL_H + +#include "matching.h" + +MATCHING_NAMESPACE_BEGIN; + +template + KDNodeStats::KDNodeStats(const Matrix &points) + : KDNode(points), pointStats_(*(new all_point_stats_type(n_points_))), changed_(true) +{ +} + +template + KDNodeStats::KDNodeStats(KDNodeStats *parent) + : KDNode(parent), pointStats_(parent->pointStats_), changed_(true) +{ +} + +template + KDNodeStats::~KDNodeStats() +{ + if (isRoot()) delete &pointStats_; +} + +template + KDNode* KDNodeStats::newNode(KDNode *parent) +{ + return new KDNodeStats((KDNodeStats*) parent); +} + +template + void KDNodeStats::setPointStats(int index, const point_stats_type &stats) +{ + pointStats_[oldIndex(index)] = stats; + ((KDNodeStats*)leaf(index))->setChanged(true); +} + +template + void KDNodeStats::setChanged(bool changed) +{ + if (changed_ == changed) return; + changed_ = changed; + if (changed && !isRoot()) + ((KDNodeStats*)parent())->setChanged(changed); +} + +template + void KDNodeStats::visit(bool init) +{ + if (!isChanged()) return; // the node is unchanged, not neccessary to proceed + if (isLeaf()) + setLeafStats(init); + else + { + for (unsigned int i = 0; i < children_.size(); i++) + { + ((KDNodeStats*) children_[i])->visit(init); + } + setNonLeafStats(init); + } + setChanged(false); // the node statistics is refreshed +} + +MATCHING_NAMESPACE_END; + +#endif // KDNODE_IMPL_H diff --git a/fastlib/trunk/contrib/tqlong/matching/matching.cc b/fastlib/trunk/contrib/tqlong/matching/matching.cc index 163faf488f..ed96af92d0 100644 --- a/fastlib/trunk/contrib/tqlong/matching/matching.cc +++ b/fastlib/trunk/contrib/tqlong/matching/matching.cc @@ -14,4 +14,19 @@ std::string toString (const Vector& v) return s.str(); } +std::string toString (const Matrix& M) +{ + std::stringstream s; + for (int i = 0; i < M.n_rows(); i++) + { + if (i > 0) s << " "; + else s << "("; + for (int j = 0; j < M.n_cols(); j++) + s << " " << M.get(i, j); + if (i < M.n_rows()-1) s << "\n"; + else s << ")"; + } + return s.str(); +} + MATCHING_NAMESPACE_END; diff --git a/fastlib/trunk/contrib/tqlong/matching/matching.h b/fastlib/trunk/contrib/tqlong/matching/matching.h index 0a870c2999..cd8149c083 100644 --- a/fastlib/trunk/contrib/tqlong/matching/matching.h +++ b/fastlib/trunk/contrib/tqlong/matching/matching.h @@ -16,6 +16,8 @@ MATCHING_NAMESPACE_BEGIN; std::string toString (const Vector& v); +std::string toString (const Matrix& v); + MATCHING_NAMESPACE_END; #endif // MATCHING_H diff --git a/fastlib/trunk/contrib/tqlong/matching/mt_test.cpp b/fastlib/trunk/contrib/tqlong/matching/mt_test.cpp index 234ac8f24a..ae697e42d4 100644 --- a/fastlib/trunk/contrib/tqlong/matching/mt_test.cpp +++ b/fastlib/trunk/contrib/tqlong/matching/mt_test.cpp @@ -7,6 +7,7 @@ #include "matching.h" #include "kdnode.h" #include "single_tree.h" +#include "auction_matching.h" //namespace po = boost::program_options; using namespace std; @@ -170,19 +171,103 @@ template <> return sqrt(s)+stats.minPrice_; } -MATCHING_NAMESPACE_END; - -int main(int argc, char** argv) +class Mat : public Matrix { - process_options(argc, argv); - - if (vm["random"].as() > 0) + std::vector price_; +public: + void Init(int rows, int cols) { - generateRandom(vm["reference"].as().c_str(), vm["query"].as().c_str()); + Matrix::Init(rows, cols); + price_ = std::vector(n_cols(), 0); + } + void setPrice(int col, double price) { price_[col] = price; } + double price(int col) const { return price_.at(col); } + double getP(int row, int col) const { return get(row, col)-price(col); } + double kBest(int row, std::vector &cols) + { + int k = (int) cols.size(); + std::vector maxs(k, -std::numeric_limits::infinity()); + for (int col = 0; col < n_cols(); col++) + { + double val = getP(row, col); + for (int o = 0; o < k; o++) if (val > maxs[o]) + { + for (int k_r = k-1; k_r > o; k_r--) + { + maxs[k_r] = maxs[k_r-1]; + cols[k_r] = cols[k_r-1]; + } + maxs[o] = val; + cols[o] = col; + break; + } + } + return 0; + } + void refresh() {} +}; + +class KDMat +{ +protected: + typedef KDNodeStats node_type; + typedef match::SingleTree SingleTree; + const Matrix &ref_, &query_; + node_type* rRoot; +public: + KDMat(const Matrix& reference, const Matrix& query) + : ref_(reference), query_(query) + { +// cout << "start 0\n"; + rRoot = new node_type(ref_); + rRoot->split(); + for (int col = 0; col < n_cols(); col++) + rRoot->setPointStats(col, 0); + rRoot->visit(true); +// cout << "done 0\n"; + } + ~KDMat() + { + delete rRoot; } - ptime time_start(second_clock::local_time()); - + int n_rows() const { return query_.n_cols(); } + int n_cols() const { return rRoot->n_points(); } + double get(int row, int col) const + { + Vector q_vec; + queryVec(row, q_vec); + return -SingleTree::distance(q_vec, *rRoot, col); + } + void setPrice(int col, PointStats price) { rRoot->setPointStats(col, price); } + PointStats price(int col) const { return rRoot->pointStats(col); } + double kBest(int row, std::vector &cols) + { +// cout << "start 1\n"; + std::vector mins((int) cols.size(), std::numeric_limits::infinity()); + Vector q_vec; + queryVec(row, q_vec); + return SingleTree::kNearestNeighbor(q_vec, *rRoot, cols, mins); +// cout << "done 1\n"; + } + void refresh() + { + rRoot->visit(false); + } + void queryVec(int row, Vector& q) const + { + query_.MakeColumnVector(row, &q); + } + void refVec(int col, Vector& r) const + { + rRoot->getPoint(col, r); + } +}; + +MATCHING_NAMESPACE_END; + +void testKDNodeStats() +{ Matrix reference, query; data::Load(vm["reference"].as().c_str(), &reference); data::Load(vm["query"].as().c_str(), &query); @@ -195,24 +280,86 @@ int main(int argc, char** argv) cout << "Done set stats\n"; rRoot->visit(true); cout << "Done visit\n"; - cout << rRoot->toString() << "\n"; +// cout << rRoot->toString() << "\n"; typedef match::SingleTree SingleTree; - SingleTree algo; + double total = 0; for (int i = 0; i < query.n_cols(); i++) { - int minIndex = -1; double minDistance = std::numeric_limits::infinity(); - Vector q, r; + Vector q; query.MakeColumnVector(i, &q); - double pruned = algo.nearestNeighbor(q, *rRoot, minIndex, minDistance); - rRoot->getPoint(minIndex, r); - cout << "q = " << match::toString(q) << " --> r = " << match::toString(r) - << " index = " << minIndex << " dist = " << minDistance - << " pruned = " << pruned << "\n"; +// int minIndex = -1; double minDistance = std::numeric_limits::infinity(); +// pruned += algo.nearestNeighbor(q, *rRoot, minIndex, minDistance); +// Vector r; +// rRoot->getPoint(minIndex, r); +// cout << "q = " << match::toString(q) << " --> r = " << match::toString(r) +// << " index = " << minIndex << " dist = " << minDistance +// << " pruned = " << pruned << "\n"; + std::vector minIndexes(2, -1); std::vector minDistances(2, std::numeric_limits::infinity()); + double pruned = SingleTree::kNearestNeighbor(q, *rRoot, minIndexes, minDistances); +// cout << "q = " << match::toString(q) << " --> " << " index = " << minIndexes[0] << " " << minIndexes[1] +// << " dist = " << minDistances[0] << " " << minDistances[1] +// << " pruned = " << pruned << "\n"; + total += pruned; } + cout << "Done finding nearest neighbors pruned = " << total << "\n"; + delete rRoot; +} + +void testMatching() +{ + int n = vm["random"].as(); + match::Mat W; + W.Init(n,n); + for (int i = 0; i < W.n_rows(); i++) + for (int j = 0; j < W.n_cols(); j++) + W.ref(i, j) = math::RandInt(0, 10); + +// cout << match::toString(W) << "\n"; + + match::AuctionMatching auction(W); + auction.doMatch(); + + for (int i = 0; i < auction.n_left(); i++) + cout << i << " --> " << auction.left(i) << " w = " << W.get(i, auction.left(i)) << "\n"; +} + +void testMatching1() +{ + Matrix reference, query; + data::Load(vm["reference"].as().c_str(), &reference); + data::Load(vm["query"].as().c_str(), &query); + + typedef match::KDNodeStats Node; + typedef match::KDMat Mat; + +// cout << match::toString(W) << "\n"; + + Mat W(reference, query); + match::AuctionMatching auction(W); + auction.doMatch(); + +// for (int i = 0; i < auction.n_left(); i++) +// cout << i << " --> " << auction.left(i) << " w = " << W.get(i, auction.left(i)) << "\n"; +} + +int main(int argc, char** argv) +{ + process_options(argc, argv); + + if (vm["random"].as() > 0) + { + generateRandom(vm["reference"].as().c_str(), vm["query"].as().c_str()); + } + + ptime time_start(second_clock::local_time()); + +// testKDNodeStats(); +// testMatching(); + testMatching1(); ptime time_end(second_clock::local_time()); time_duration duration(time_end - time_start); diff --git a/fastlib/trunk/contrib/tqlong/matching/single_tree.h b/fastlib/trunk/contrib/tqlong/matching/single_tree.h index 23f632192b..59e4a4cc4a 100644 --- a/fastlib/trunk/contrib/tqlong/matching/single_tree.h +++ b/fastlib/trunk/contrib/tqlong/matching/single_tree.h @@ -15,13 +15,14 @@ public: typedef P point_type; typedef T node_type; - double nearestNeighbor(const point_type& q, node_type& ref, int &minIndex, double &minDistance); + static double nearestNeighbor(const point_type& q, node_type& ref, int &minIndex, double &minDistance); + static double kNearestNeighbor(const point_type& q, node_type& ref, std::vector& minIndexes, std::vector& minDistances); // distance of q from point index in ref - double distance(const point_type& q, node_type& ref, int index); + static double distance(const point_type& q, node_type& ref, int index); // distance of q to bounding box of ref - double distance(const point_type& q, node_type& ref); + static double distance(const point_type& q, node_type& ref); }; template @@ -41,7 +42,7 @@ template if (val < minDistance) { minDistance = val; - minIndex = i; + minIndex = ref.index(i); // return index from root view } } return 0; @@ -69,6 +70,49 @@ template // } } +template + double SingleTree::kNearestNeighbor(const point_type &q, node_type &ref, + std::vector &minIndexes, std::vector &minDistances) +{ + int k = (int) minIndexes.size(); + if (distance(q, ref) >= minDistances[k-1]) + { + return ref.n_points(); + } + if (ref.isLeaf()) + { + for (int i = 0; i < ref.n_points(); i++) + { + double val = distance(q, ref, i); + for (int order = 0; order < k; order++) + { + if (val < minDistances[order]) + { + for (int k_reverse = k-1; k_reverse > order; k_reverse--) + { + minDistances[k_reverse] = minDistances[k_reverse-1]; + minIndexes[k_reverse] = minIndexes[k_reverse-1]; + } + minDistances[order] = val; + minIndexes[order] = ref.index(i); // return index from root view + break; + } + } + } + return 0; + } + else + { + double pruned = 0; + for (int i = 0; i < ref.n_children(); i++) + { + node_type* child = (node_type*)ref.child(i); + pruned += kNearestNeighbor(q, *child, minIndexes, minDistances); + } + return pruned; + } +} + MATCHING_NAMESPACE_END; #endif // SINGLE_TREE_H