matching: auction templates

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

fastlib/trunk/contrib/tqlong/matching/auction_matching.h

@@ -0,0 +1,120 @@
+#ifndef AUCTION_MATCHING_H
+#define AUCTION_MATCHING_H
+
+#include <algorithm>
+#include "matching.h"
+
+MATCHING_NAMESPACE_BEGIN;
+
+template <typename W>
+    class AuctionMatching
+{
+public:
+  typedef W                weight_type;
+protected:
+  weight_type&        weight_;
+  double              epsilon_;
+  std::vector<int>    left_, right_, bidders_;
+  std::vector<double> 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 <typename W>
+    AuctionMatching<W>::AuctionMatching(weight_type &weight)
+      : weight_(weight), left_(n_left()), right_(n_right()),
+        bidders_(n_right()), bids_(n_right())
+{
+}
+
+template <typename W>
+    void AuctionMatching<W>::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 <typename W>
+    void AuctionMatching<W>::clearMatches()
+{
+  std::fill(left_.begin(), left_.end(), -1);
+  std::fill(right_.begin(), right_.end(), -1);
+}
+
+template <typename W>
+    void AuctionMatching<W>::forwardAuction(double &pruned, double &total)
+{
+  pruned = total = 0;
+  while (1)
+  {
+    bool doneMatching = true;
+    std::fill(bids_.begin(), bids_.end(), -std::numeric_limits<double>::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<int> 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 <typename W>
+    void AuctionMatching<W>::placeBid(int l, int r, double price)
+{
+  if (bids_[r] < price)
+  {
+    bids_[r] = price;
+    bidders_[r] = l;
+  }
+}
+
+template <typename W>
+    void AuctionMatching<W>::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

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);

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 <typename P, typename N>
-    KDNodeStats<P,N>::KDNodeStats(const Matrix &points)
-      : KDNode(points), pointStats_(*(new all_point_stats_type(n_points_))), changed_(true)
-{
-}
-
-template <typename P, typename N>
-    KDNodeStats<P,N>::KDNodeStats(KDNodeStats *parent)
-      : KDNode(parent), pointStats_(parent->pointStats_), changed_(true)
-{
-}
-
-template <typename P, typename N>
-    KDNodeStats<P,N>::~KDNodeStats()
-{
-  if (isRoot()) delete &pointStats_;
-}
-
-template <typename P, typename N>
-    KDNode* KDNodeStats<P,N>::newNode(KDNode *parent)
-{
-  return new KDNodeStats((KDNodeStats*) parent);
-}
-
-template <typename P, typename N>
-    void KDNodeStats<P,N>::setPointStats(int index, const point_stats_type &stats)
-{
-  pointStats_[oldIndex(index)] = stats;
-  ((KDNodeStats*)leaf(index))->setChanged(true);
-}
-
-template <typename P, typename N>
-    void KDNodeStats<P,N>::setChanged(bool changed)
-{
-  if (changed_ == changed) return;
-  changed_ = changed;
-  if (changed && !isRoot())
-    ((KDNodeStats*)parent())->setChanged(changed);
-}
-
-template <typename P, typename N>
-    void KDNodeStats<P,N>::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

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 <typename P, typename N>
+    KDNodeStats<P,N>::KDNodeStats(const Matrix &points)
+      : KDNode(points), pointStats_(*(new all_point_stats_type(n_points_))), changed_(true)
+{
+}
+
+template <typename P, typename N>
+    KDNodeStats<P,N>::KDNodeStats(KDNodeStats *parent)
+      : KDNode(parent), pointStats_(parent->pointStats_), changed_(true)
+{
+}
+
+template <typename P, typename N>
+    KDNodeStats<P,N>::~KDNodeStats()
+{
+  if (isRoot()) delete &pointStats_;
+}
+
+template <typename P, typename N>
+    KDNode* KDNodeStats<P,N>::newNode(KDNode *parent)
+{
+  return new KDNodeStats((KDNodeStats*) parent);
+}
+
+template <typename P, typename N>
+    void KDNodeStats<P,N>::setPointStats(int index, const point_stats_type &stats)
+{
+  pointStats_[oldIndex(index)] = stats;
+  ((KDNodeStats*)leaf(index))->setChanged(true);
+}
+
+template <typename P, typename N>
+    void KDNodeStats<P,N>::setChanged(bool changed)
+{
+  if (changed_ == changed) return;
+  changed_ = changed;
+  if (changed && !isRoot())
+    ((KDNodeStats*)parent())->setChanged(changed);
+}
+
+template <typename P, typename N>
+    void KDNodeStats<P,N>::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

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;

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

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<int>() > 0)
+  std::vector<double> price_;
+public:
+  void Init(int rows, int cols)
   {
-    generateRandom(vm["reference"].as<string>().c_str(), vm["query"].as<string>().c_str());
+    Matrix::Init(rows, cols);
+    price_ = std::vector<double>(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<int> &cols)
+  {
+    int k = (int) cols.size();
+    std::vector<double> maxs(k, -std::numeric_limits<double>::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<PointStats, NodeStats>    node_type;
+  typedef match::SingleTree<Vector, node_type>  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<int> &cols)
+  {
+//    cout << "start 1\n";
+    std::vector<double> mins((int) cols.size(), std::numeric_limits<double>::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<string>().c_str(), &reference);
   data::Load(vm["query"].as<string>().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<Vector, Node> SingleTree;
-  SingleTree algo;
 
+  double total = 0;
   for (int i = 0; i < query.n_cols(); i++)
   {
-    int minIndex = -1; double minDistance = std::numeric_limits<double>::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<double>::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<int> minIndexes(2, -1); std::vector<double> minDistances(2, std::numeric_limits<double>::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<int>();
+  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<match::Mat> 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<string>().c_str(), &reference);
+  data::Load(vm["query"].as<string>().c_str(), &query);
+
+  typedef match::KDNodeStats<match::PointStats, match::NodeStats> Node;
+  typedef match::KDMat Mat;
+
+//  cout << match::toString(W) << "\n";
+
+  Mat W(reference, query);
+  match::AuctionMatching<Mat> 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<int>() > 0)
+  {
+    generateRandom(vm["reference"].as<string>().c_str(), vm["query"].as<string>().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);

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<int>& minIndexes, std::vector<double>& 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 <typename P, typename T>
@@ -41,7 +42,7 @@ template <typename P, typename T>
       if (val < minDistance)
       {
         minDistance = val;
-        minIndex = i;
+        minIndex = ref.index(i); // return index from root view
       }
     }
     return 0;
@@ -69,6 +70,49 @@ template <typename P, typename T>
 //  }
 }
 
+template <typename P, typename T>
+    double SingleTree<P,T>::kNearestNeighbor(const point_type &q, node_type &ref,
+                                             std::vector<int> &minIndexes, std::vector<double> &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