From 15e84d930a495069eecf76a1d155dc913f06dd97 Mon Sep 17 00:00:00 2001
From: Patrick Brosi <brosi@informatik.uni-freiburg.de>
Date: Fri, 21 Jan 2022 23:58:49 +0100
Subject: [PATCH] update segmentation in shapevl, add --unique mode

---
 src/pfaedle/router/Router.h                   | 17 ++--
 src/pfaedle/router/Router.tpp                 | 15 ++--
 src/pfaedle/router/RoutingAttrs.h             |  9 ++-
 src/pfaedle/router/ShapeBuilder.cpp           | 12 +--
 src/pfaedle/router/ShapeBuilder.h             | 14 ++--
 src/pfaedle/router/TripTrie.h                 | 14 ++--
 .../router/{TripTrie.cpp => TripTrie.tpp}     | 30 +++++---
 src/shapevl/Collector.cpp                     | 75 +++++++++++-------
 src/shapevl/Collector.h                       |  1 +
 src/shapevl/ShapevlMain.cpp                   | 77 ++++++++++++++-----
 src/util/geo/Geo.h                            |  8 +-
 11 files changed, 177 insertions(+), 95 deletions(-)
 rename src/pfaedle/router/{TripTrie.cpp => TripTrie.tpp} (88%)

diff --git a/src/pfaedle/router/Router.h b/src/pfaedle/router/Router.h
index 4ec8cb7..630d456 100644
--- a/src/pfaedle/router/Router.h
+++ b/src/pfaedle/router/Router.h
@@ -6,11 +6,11 @@
 #define PFAEDLE_ROUTER_ROUTER_H_
 
 #include <limits>
+#include <map>
 #include <set>
 #include <stack>
 #include <string>
 #include <unordered_map>
-#include <map>
 #include <utility>
 #include <vector>
 #include "pfaedle/Def.h"
@@ -52,12 +52,10 @@ typedef std::vector<std::pair<std::pair<size_t, size_t>, uint32_t>> CostMatrix;
 class Router {
  public:
   virtual ~Router() = default;
-  virtual std::map<size_t, EdgeListHops> route(const TripTrie* trie,
-                                               const EdgeCandMap& ecm,
-                                               const RoutingOpts& rOpts,
-                                               const osm::Restrictor& rest,
-                                               HopCache* hopCache,
-                                               bool noFastHops) const = 0;
+  virtual std::map<size_t, EdgeListHops> route(
+      const TripTrie<pfaedle::gtfs::Trip>* trie, const EdgeCandMap& ecm,
+      const RoutingOpts& rOpts, const osm::Restrictor& rest, HopCache* hopCache,
+      bool noFastHops) const = 0;
 };
 
 /*
@@ -69,8 +67,9 @@ class RouterImpl : public Router {
  public:
   // Find the most likely path through the graph for a trip trie.
   virtual std::map<size_t, EdgeListHops> route(
-      const TripTrie* trie, const EdgeCandMap& ecm, const RoutingOpts& rOpts,
-      const osm::Restrictor& rest, HopCache* hopCache, bool noFastHops) const;
+      const TripTrie<pfaedle::gtfs::Trip>* trie, const EdgeCandMap& ecm,
+      const RoutingOpts& rOpts, const osm::Restrictor& rest, HopCache* hopCache,
+      bool noFastHops) const;
 
  private:
   void hops(const EdgeCandGroup& from, const EdgeCandGroup& to,
diff --git a/src/pfaedle/router/Router.tpp b/src/pfaedle/router/Router.tpp
index 9b35f51..a62d11a 100644
--- a/src/pfaedle/router/Router.tpp
+++ b/src/pfaedle/router/Router.tpp
@@ -9,21 +9,22 @@
 #define omp_get_num_procs() 1
 #endif
 
-#include <unordered_map>
-#include <map>
-#include <vector>
-#include <utility>
-#include <set>
 #include <limits>
+#include <map>
+#include <set>
 #include <stack>
+#include <unordered_map>
+#include <utility>
+#include <vector>
 
 using util::graph::EDijkstra;
 
 // _____________________________________________________________________________
 template <typename TW>
 std::map<size_t, EdgeListHops> RouterImpl<TW>::route(
-    const TripTrie* trie, const EdgeCandMap& ecm, const RoutingOpts& rOpts,
-    const osm::Restrictor& rest, HopCache* hopCache, bool noFastHops) const {
+    const TripTrie<pfaedle::gtfs::Trip>* trie, const EdgeCandMap& ecm,
+    const RoutingOpts& rOpts, const osm::Restrictor& rest, HopCache* hopCache,
+    bool noFastHops) const {
   std::map<size_t, EdgeListHops> ret;
 
   // the current node costs in our DAG
diff --git a/src/pfaedle/router/RoutingAttrs.h b/src/pfaedle/router/RoutingAttrs.h
index c1b001f..3f7965b 100644
--- a/src/pfaedle/router/RoutingAttrs.h
+++ b/src/pfaedle/router/RoutingAttrs.h
@@ -5,9 +5,9 @@
 #ifndef PFAEDLE_ROUTER_ROUTINGATTRS_H_
 #define PFAEDLE_ROUTER_ROUTINGATTRS_H_
 
+#include <string>
 #include <unordered_map>
 #include <vector>
-#include <string>
 #include "pfaedle/statsimi-classifier/StatsimiClassifier.h"
 #include "pfaedle/trgraph/EdgePL.h"
 
@@ -30,6 +30,13 @@ inline bool operator<(const LineSimilarity& a, const LineSimilarity& b) {
 struct RoutingAttrs {
   RoutingAttrs()
       : lineFrom(""), lineTo(), shortName(""), classifier(0), _simiCache() {}
+  RoutingAttrs(const std::string& shortName, const std::string& lineFrom,
+               const std::string& lineTo)
+      : lineFrom(lineFrom),
+        lineTo({lineTo}),
+        shortName(shortName),
+        classifier(0),
+        _simiCache() {}
   std::string lineFrom;
   std::vector<std::string> lineTo;
   std::string shortName;
diff --git a/src/pfaedle/router/ShapeBuilder.cpp b/src/pfaedle/router/ShapeBuilder.cpp
index 3615ed1..1f69696 100644
--- a/src/pfaedle/router/ShapeBuilder.cpp
+++ b/src/pfaedle/router/ShapeBuilder.cpp
@@ -328,14 +328,15 @@ std::pair<std::vector<LINE>, Stats> ShapeBuilder::shapeL(Trip* trip) {
 
 // _____________________________________________________________________________
 std::map<size_t, pfaedle::router::EdgeListHops> ShapeBuilder::route(
-    const TripTrie* trie, const EdgeCandMap& ecm, HopCache* hopCache) const {
+    const TripTrie<pfaedle::gtfs::Trip>* trie, const EdgeCandMap& ecm,
+    HopCache* hopCache) const {
   return _router->route(trie, ecm, _motCfg.routingOpts, *_restr, hopCache,
                         _cfg.noFastHops);
 }
 
 // _____________________________________________________________________________
 std::map<size_t, EdgeListHops> ShapeBuilder::shapeify(
-    const TripTrie* trie, HopCache* hopCache) const {
+    const TripTrie<pfaedle::gtfs::Trip>* trie, HopCache* hopCache) const {
   LOG(VDEBUG) << "Map-matching trie " << trie;
 
   // TODO(patrick): assumes the trie is not empty, check this!
@@ -362,7 +363,7 @@ EdgeListHops ShapeBuilder::shapeify(Trip* trip) {
               << trip->getRoute()->getType() << "(sn=" << trip->getShortname()
               << ", rsn=" << trip->getRoute()->getShortName()
               << ", rln=" << trip->getRoute()->getLongName() << ")";
-  TripTrie trie;
+  TripTrie<pfaedle::gtfs::Trip> trie;
   trie.addTrip(trip, getRAttrs(trip),
                _motCfg.routingOpts.transPenMethod == "timenorm", false);
   const auto& routes = route(&trie, getECM(&trie), 0);
@@ -748,7 +749,8 @@ std::vector<double> ShapeBuilder::getTransDists(Trip* trip) const {
 }
 
 // _____________________________________________________________________________
-EdgeCandMap ShapeBuilder::getECM(const TripTrie* trie) const {
+EdgeCandMap ShapeBuilder::getECM(
+    const TripTrie<pfaedle::gtfs::Trip>* trie) const {
   EdgeCandMap ecm(trie->getNds().size());
 
   for (size_t nid = 1; nid < trie->getNds().size(); nid++) {
@@ -1146,7 +1148,7 @@ void ShapeBuilder::shapeWorker(
     if (!_cfg.noHopCache) hopCache = &hopCacheLoc;
 
     for (size_t i = 0; i < forest.size(); i++) {
-      const TripTrie* trie = &(forest[i]);
+      const TripTrie<pfaedle::gtfs::Trip>* trie = &(forest[i]);
       const auto& hops = shapeify(trie, hopCache);
 
       for (const auto& leaf : trie->getNdTrips()) {
diff --git a/src/pfaedle/router/ShapeBuilder.h b/src/pfaedle/router/ShapeBuilder.h
index 3e243ad..662ca59 100644
--- a/src/pfaedle/router/ShapeBuilder.h
+++ b/src/pfaedle/router/ShapeBuilder.h
@@ -30,7 +30,7 @@
 namespace pfaedle {
 namespace router {
 
-typedef std::vector<TripTrie> TripForest;
+typedef std::vector<TripTrie<pfaedle::gtfs::Trip>> TripForest;
 typedef std::map<router::RoutingAttrs, TripForest> TripForests;
 typedef std::pair<const ad::cppgtfs::gtfs::Stop*,
                   const ad::cppgtfs::gtfs::Stop*>
@@ -64,8 +64,8 @@ class ShapeBuilder {
   // shape single trip
   std::pair<std::vector<LINE>, Stats> shapeL(pfaedle::gtfs::Trip* trip);
 
-  std::map<size_t, EdgeListHops> shapeify(const TripTrie* trie,
-                                          HopCache* hopCache) const;
+  std::map<size_t, EdgeListHops> shapeify(
+      const TripTrie<pfaedle::gtfs::Trip>* trie, HopCache* hopCache) const;
   EdgeListHops shapeify(pfaedle::gtfs::Trip* trip);
 
   const trgraph::Graph* getGraph() const;
@@ -112,14 +112,14 @@ class ShapeBuilder {
 
   EdgeCandGroup getEdgCands(const ad::cppgtfs::gtfs::Stop* s) const;
 
-  router::EdgeCandMap getECM(const TripTrie* trie) const;
+  router::EdgeCandMap getECM(const TripTrie<pfaedle::gtfs::Trip>* trie) const;
   std::vector<double> getTransTimes(pfaedle::gtfs::Trip* trip) const;
   std::vector<double> getTransDists(pfaedle::gtfs::Trip* trip) const;
   const router::RoutingAttrs& getRAttrs(const pfaedle::gtfs::Trip* trip) const;
   const router::RoutingAttrs& getRAttrs(const pfaedle::gtfs::Trip* trip);
-  std::map<size_t, router::EdgeListHops> route(const TripTrie* trie,
-                                               const EdgeCandMap& ecm,
-                                               HopCache* hopCache) const;
+  std::map<size_t, router::EdgeListHops> route(
+      const TripTrie<pfaedle::gtfs::Trip>* trie, const EdgeCandMap& ecm,
+      HopCache* hopCache) const;
   double emWeight(double mDist) const;
 
   void buildCandCache(const TripForests& clusters);
diff --git a/src/pfaedle/router/TripTrie.h b/src/pfaedle/router/TripTrie.h
index e03217f..dfcfde5 100644
--- a/src/pfaedle/router/TripTrie.h
+++ b/src/pfaedle/router/TripTrie.h
@@ -31,26 +31,27 @@ struct TripTrieNd {
   RoutingAttrs rAttrs;
 };
 
+template <typename TRIP>
 class TripTrie {
  public:
   // init node 0, this is the first decision node
   TripTrie() : _nds(1) {}
-  bool addTrip(pfaedle::gtfs::Trip* trip, const RoutingAttrs& rAttrs,
+  bool addTrip(TRIP* trip, const RoutingAttrs& rAttrs,
                bool timeEx, bool degen);
 
   const std::vector<TripTrieNd>& getNds() const;
   const TripTrieNd& getNd(size_t nid) const;
 
   void toDot(std::ostream& os, const std::string& rootName, size_t gid) const;
-  const std::map<size_t, std::vector<pfaedle::gtfs::Trip*>>& getNdTrips() const;
+  const std::map<size_t, std::vector<TRIP*>>& getNdTrips() const;
 
  private:
   std::vector<TripTrieNd> _nds;
-  std::map<pfaedle::gtfs::Trip*, size_t> _tripNds;
-  std::map<size_t, std::vector<pfaedle::gtfs::Trip*>> _ndTrips;
+  std::map<TRIP*, size_t> _tripNds;
+  std::map<size_t, std::vector<TRIP*>> _ndTrips;
 
-  bool add(pfaedle::gtfs::Trip* trip, const RoutingAttrs& rAttrs, bool timeEx);
-  size_t get(pfaedle::gtfs::Trip* trip, bool timeEx);
+  bool add(TRIP* trip, const RoutingAttrs& rAttrs, bool timeEx);
+  size_t get(TRIP* trip, bool timeEx);
 
   size_t getMatchChild(size_t parentNid, const std::string& stopName,
                        const std::string& platform, POINT pos, int time,
@@ -59,6 +60,7 @@ class TripTrie {
                 const POINT& pos, int time, bool arr, size_t parent);
 };
 
+#include "pfaedle/router/TripTrie.tpp"
 }  // namespace router
 }  // namespace pfaedle
 
diff --git a/src/pfaedle/router/TripTrie.cpp b/src/pfaedle/router/TripTrie.tpp
similarity index 88%
rename from src/pfaedle/router/TripTrie.cpp
rename to src/pfaedle/router/TripTrie.tpp
index 470ee09..4b2b87e 100644
--- a/src/pfaedle/router/TripTrie.cpp
+++ b/src/pfaedle/router/TripTrie.tpp
@@ -9,13 +9,13 @@
 #include "ad/cppgtfs/gtfs/Feed.h"
 #include "pfaedle/gtfs/Feed.h"
 #include "pfaedle/gtfs/StopTime.h"
-#include "pfaedle/router/TripTrie.h"
 
 using pfaedle::gtfs::Trip;
 using pfaedle::router::TripTrie;
 
 // _____________________________________________________________________________
-bool TripTrie::addTrip(pfaedle::gtfs::Trip* trip, const RoutingAttrs& rAttrs,
+template <typename TRIP>
+bool TripTrie<TRIP>::addTrip(TRIP* trip, const RoutingAttrs& rAttrs,
                        bool timeEx, bool degen) {
   if (!degen) return add(trip, rAttrs, timeEx);
 
@@ -31,7 +31,8 @@ bool TripTrie::addTrip(pfaedle::gtfs::Trip* trip, const RoutingAttrs& rAttrs,
 }
 
 // _____________________________________________________________________________
-bool TripTrie::add(pfaedle::gtfs::Trip* trip, const RoutingAttrs& rAttrs,
+template <typename TRIP>
+bool TripTrie<TRIP>::add(TRIP* trip, const RoutingAttrs& rAttrs,
                    bool timeEx) {
   if (trip->getStopTimes().size() == 0) return false;
 
@@ -92,7 +93,8 @@ bool TripTrie::add(pfaedle::gtfs::Trip* trip, const RoutingAttrs& rAttrs,
 }
 
 // _____________________________________________________________________________
-size_t TripTrie::get(pfaedle::gtfs::Trip* trip, bool timeEx) {
+template <typename TRIP>
+size_t TripTrie<TRIP>::get(TRIP* trip, bool timeEx) {
   if (trip->getStopTimes().size() == 0) return false;
 
   int startSecs = trip->getStopTimes().front().getDepartureTime().seconds();
@@ -137,7 +139,8 @@ size_t TripTrie::get(pfaedle::gtfs::Trip* trip, bool timeEx) {
 }
 
 // _____________________________________________________________________________
-size_t TripTrie::insert(const ad::cppgtfs::gtfs::Stop* stop,
+template <typename TRIP>
+size_t TripTrie<TRIP>::insert(const ad::cppgtfs::gtfs::Stop* stop,
                         const RoutingAttrs& rAttrs, const POINT& pos, int time,
                         bool arr, size_t parent) {
   _nds.emplace_back(TripTrieNd{stop,
@@ -158,12 +161,14 @@ size_t TripTrie::insert(const ad::cppgtfs::gtfs::Stop* stop,
 }
 
 // _____________________________________________________________________________
-const std::vector<pfaedle::router::TripTrieNd>& TripTrie::getNds() const {
+template <typename TRIP>
+const std::vector<pfaedle::router::TripTrieNd>& TripTrie<TRIP>::getNds() const {
   return _nds;
 }
 
 // _____________________________________________________________________________
-size_t TripTrie::getMatchChild(size_t parentNid, const std::string& stopName,
+template <typename TRIP>
+size_t TripTrie<TRIP>::getMatchChild(size_t parentNid, const std::string& stopName,
                                const std::string& platform, POINT pos, int time,
                                bool timeEx) const {
   for (size_t child : _nds[parentNid].childs) {
@@ -179,7 +184,8 @@ size_t TripTrie::getMatchChild(size_t parentNid, const std::string& stopName,
 }
 
 // _____________________________________________________________________________
-void TripTrie::toDot(std::ostream& os, const std::string& rootName,
+template <typename TRIP>
+void TripTrie<TRIP>::toDot(std::ostream& os, const std::string& rootName,
                      size_t gid) const {
   os << "digraph triptrie" << gid << " {";
 
@@ -208,12 +214,14 @@ void TripTrie::toDot(std::ostream& os, const std::string& rootName,
 }
 
 // _____________________________________________________________________________
-const std::map<size_t, std::vector<pfaedle::gtfs::Trip*>>&
-TripTrie::getNdTrips() const {
+template <typename TRIP>
+const std::map<size_t, std::vector<TRIP*>>&
+TripTrie<TRIP>::getNdTrips() const {
   return _ndTrips;
 }
 
 // _____________________________________________________________________________
-const pfaedle::router::TripTrieNd& TripTrie::getNd(size_t nid) const {
+template <typename TRIP>
+const pfaedle::router::TripTrieNd& TripTrie<TRIP>::getNd(size_t nid) const {
   return _nds[nid];
 }
diff --git a/src/shapevl/Collector.cpp b/src/shapevl/Collector.cpp
index 4d6f98c..c6c78c9 100644
--- a/src/shapevl/Collector.cpp
+++ b/src/shapevl/Collector.cpp
@@ -27,6 +27,7 @@ using util::geo::output::GeoJsonOutput;
 double Collector::add(const Trip* oldT, const Shape* oldS, const Trip* newT,
                       const Shape* newS) {
   // This adds a new trip with a new shape to our evaluation.
+  // if (oldT->getId() != "pse-a779ac00") return 0;
 
   _trips++;
 
@@ -36,7 +37,6 @@ double Collector::add(const Trip* oldT, const Shape* oldS, const Trip* newT,
     return 0;
   }
 
-
   for (auto st : oldT->getStopTimes()) {
     if (st.getShapeDistanceTravelled() < 0) {
       // we cannot safely compare trips without shape dist travelled
@@ -82,7 +82,6 @@ double Collector::add(const Trip* oldT, const Shape* oldS, const Trip* newT,
     return 0;
   }
 
-
   std::vector<std::pair<double, double>> newLenDists;
   std::vector<std::pair<double, double>> oldLenDists;
 
@@ -111,8 +110,10 @@ double Collector::add(const Trip* oldT, const Shape* oldS, const Trip* newT,
   double f = util::geo::webMercDistFactor(oldLCut.front());
 
   // roughly half a meter
-  auto oldLCutS = util::geo::simplify(oldLCut, f * (0.5 / util::geo::M_PER_DEG));
-  auto newLCutS = util::geo::simplify(newLCut, f * (0.5 / util::geo::M_PER_DEG));
+  auto oldLCutS =
+      util::geo::simplify(oldLCut, f * (0.5 / util::geo::M_PER_DEG));
+  auto newLCutS =
+      util::geo::simplify(newLCut, f * (0.5 / util::geo::M_PER_DEG));
 
   auto old = _dCache.find(oldLCutS);
   if (old != _dCache.end()) {
@@ -156,6 +157,7 @@ double Collector::add(const Trip* oldT, const Shape* oldS, const Trip* newT,
   if (AN <= 0.0001) _an0++;
   if (AN <= 0.05) _an5++;
   if (AN <= 0.1) _an10++;
+  if (AN <= 0.2) _an20++;
   if (AN <= 0.3) _an30++;
   if (AN <= 0.5) _an50++;
   if (AN <= 0.7) _an70++;
@@ -167,9 +169,17 @@ double Collector::add(const Trip* oldT, const Shape* oldS, const Trip* newT,
               << totL << " = " << AL << " d_f = " << avgFd;
 
   if (_reportOut) {
+    (*_reportOut) << std::fixed << std::setprecision(6);
     (*_reportOut) << oldT->getId() << "\t" << AN << "\t" << AL << "\t" << avgFd
                   << "\t" << util::geo::getWKT(oldSegs) << "\t"
-                  << util::geo::getWKT(newSegs) << "\n";
+                  << util::geo::getWKT(newSegs) << "\t" << oldT->getRoute()->getShortName() << "\t";
+
+    for (const auto& st : oldT->getStopTimes()) {
+      (*_reportOut) << st.getStop()->getName() << "\t"
+                    << st.getStop()->getLat() << "\t"
+                    << st.getStop()->getLng() << "\t";
+    }
+(*_reportOut) << "\n";
   }
 
   return avgFd;
@@ -189,33 +199,32 @@ std::vector<LINE> Collector::segmentize(
   size_t i = 0;
   for (auto st : t->getStopTimes()) {
     cuts.push_back(std::pair<POINT, double>(
-          {st.getStop()->getLng(),
-                                              st.getStop()->getLat()},
+        {st.getStop()->getLng(), st.getStop()->getLat()},
         st.getShapeDistanceTravelled()));
     i++;
   }
 
-  // get first half of geometry, and search for start point there!
-  size_t before = std::upper_bound(dists.begin(), dists.end(), cuts[1].second) -
+
+  size_t to = std::upper_bound(dists.begin(), dists.end(), cuts[0].second) -
                   dists.begin();
-  if (before + 1 > shape.size()) before = shape.size() - 1;
-  assert(shape.begin() + before + 1 <= shape.end());
-  POLYLINE l(LINE(shape.begin(), shape.begin() + before + 1));
-  auto lastLp = l.projectOn(cuts.front().first);
+  if (to >= dists.size()) to = dists.size() - 1;
+  double progr = (cuts[0].second - dists[to - 1]) / (dists[to] - dists[to - 1]);
+  auto lastP = shape[to - 1];
+  lastP.setX(lastP.getX() + progr * util::geo::dist(shape[to-1], shape[to]));
+  lastP.setY(lastP.getY() + progr * util::geo::dist(shape[to-1], shape[to]));
 
   for (size_t i = 1; i < cuts.size(); i++) {
-    size_t before = shape.size();
-    if (i < cuts.size() - 1 && cuts[i + 1].second > -0.5) {
-      before =
-          std::upper_bound(dists.begin(), dists.end(), cuts[i + 1].second) -
-          dists.begin();
-    }
 
-    POLYLINE afterPl(LINE(shape.begin(), shape.begin() + before));
+    size_t to = std::upper_bound(dists.begin(), dists.end(), cuts[i].second) -
+                    dists.begin();
+  if (to >= dists.size()) to = dists.size() - 1;
+    double progr = (cuts[i].second - dists[to - 1]) / (dists[to] - dists[to - 1]);
+    // std::cout << t->getId() << ": " << dists[to] << " vs " << cuts[i].second << ", " << dists[to - 1] << " (" << progr << ")" << std::endl;
+    auto curP = shape[to - 1];
+    curP.setX(curP.getX() + progr * util::geo::dist(shape[to-1], shape[to]));
+    curP.setY(curP.getY() + progr * util::geo::dist(shape[to-1], shape[to]));
 
-    auto curLp = afterPl.projectOnAfter(cuts[i].first, lastLp.lastIndex);
-
-    auto curL = pl.getSegment(lastLp, curLp).getLine();
+    auto curL = pl.getSegment(lastP, curP).getLine();
 
     double dist =
         util::geo::haversine(t->getStopTimes()[i - 1].getStop()->getLat(),
@@ -226,7 +235,7 @@ std::vector<LINE> Collector::segmentize(
     lenDist.push_back({dist, len});
 
     ret.push_back(curL);
-    lastLp = curLp;
+    lastP = curP;
   }
 
   return ret;
@@ -275,6 +284,10 @@ void Collector::printShortStats(std::ostream* os) const {
               static_cast<double>(_results.size())) *
                  100
           << ",";
+    (*os) << (static_cast<double>(_an20) /
+              static_cast<double>(_results.size())) *
+                 100
+          << ",";
     (*os) << (static_cast<double>(_an30) /
               static_cast<double>(_results.size())) *
                  100
@@ -332,6 +345,12 @@ void Collector::printStats(std::ostream* os) const {
                  100
           << " %"
           << "\n";
+    (*os) << "  an-20: "
+          << (static_cast<double>(_an20) /
+              static_cast<double>(_results.size())) *
+                 100
+          << " %"
+          << "\n";
     (*os) << "  acc-30: "
           << (static_cast<double>(_an30) /
               static_cast<double>(_results.size())) *
@@ -398,6 +417,8 @@ std::map<string, double> Collector::getStats() {
       (static_cast<double>(_an5) / static_cast<double>(_results.size())) * 100;
   stats["an-10"] =
       (static_cast<double>(_an10) / static_cast<double>(_results.size())) * 100;
+  stats["an-20"] =
+      (static_cast<double>(_an20) / static_cast<double>(_results.size())) * 100;
   stats["an-30"] =
       (static_cast<double>(_an30) / static_cast<double>(_results.size())) * 100;
   stats["an-50"] =
@@ -419,7 +440,7 @@ std::pair<size_t, double> Collector::getDa(const std::vector<LINE>& a,
   // convert (roughly) to degrees
   double SEGL = 25 / util::geo::M_PER_DEG;
 
-  double MAX = 50;
+  double MAX = 100;
 
   for (size_t i = 0; i < a.size(); i++) {
     double fdMeter = 0;
@@ -440,8 +461,8 @@ std::pair<size_t, double> Collector::getDa(const std::vector<LINE>& a,
         _dACache[aSimpl][bSimpl] = fdMeter;
       }
     } else {
-        fdMeter = util::geo::frechetDistHav(aSimpl, bSimpl, SEGL);
-        _dACache[aSimpl][bSimpl] = fdMeter;
+      fdMeter = util::geo::frechetDistHav(aSimpl, bSimpl, SEGL);
+      _dACache[aSimpl][bSimpl] = fdMeter;
     }
 
     if (fdMeter >= MAX) {
diff --git a/src/shapevl/Collector.h b/src/shapevl/Collector.h
index ad285e9..724f00d 100644
--- a/src/shapevl/Collector.h
+++ b/src/shapevl/Collector.h
@@ -107,6 +107,7 @@ class Collector {
   size_t _an0;
   size_t _an5;
   size_t _an10;
+  size_t _an20;
   size_t _an30;
   size_t _an50;
   size_t _an70;
diff --git a/src/shapevl/ShapevlMain.cpp b/src/shapevl/ShapevlMain.cpp
index 81f1f77..9e8c642 100644
--- a/src/shapevl/ShapevlMain.cpp
+++ b/src/shapevl/ShapevlMain.cpp
@@ -10,11 +10,14 @@
 #include <thread>
 #include <vector>
 #include "ad/cppgtfs/Parser.h"
+#include "pfaedle/router/TripTrie.h"
 #include "shapevl/Collector.h"
 #include "util/Misc.h"
 #include "util/json/Writer.h"
 #include "util/log/Log.h"
 
+using pfaedle::router::TripTrie;
+
 std::atomic<int> count(0);
 
 // _____________________________________________________________________________
@@ -37,7 +40,7 @@ void printHelp(int argc, char** argv) {
 void eval(const std::vector<std::string>* paths,
           std::vector<pfaedle::eval::Collector>* colls,
           const std::set<Route::TYPE>* mots,
-          const ad::cppgtfs::gtfs::Feed* evalFeed) {
+          const ad::cppgtfs::gtfs::Feed* evalFeed, bool unique) {
   while (1) {
     int myFeed = count-- - 1;
     if (myFeed < 0) return;
@@ -54,18 +57,59 @@ void eval(const std::vector<std::string>* paths,
       exit(1);
     }
 
+    std::vector<ad::cppgtfs::gtfs::Trip*> trips;
+
+    if (unique) {
+      std::map<const ad::cppgtfs::gtfs::Route*,
+               std::vector<TripTrie<ad::cppgtfs::gtfs::Trip>>>
+          forest;
+      for (auto t : evalFeed->getTrips()) {
+        auto& subForest = forest[t.second->getRoute()];
+        bool ins = false;
+        for (auto& trie : subForest) {
+          if (trie.addTrip(t.second,
+                           pfaedle::router::RoutingAttrs{
+                               t.second->getRoute()->getId(), "", ""},
+                           false, false)) {
+            ins = true;
+            break;
+          }
+        }
+
+        if (!ins) {
+          subForest.resize(subForest.size() + 1);
+          subForest.back().addTrip(t.second,
+                                   pfaedle::router::RoutingAttrs{
+                                       t.second->getRoute()->getId(), "", ""},
+                                   false, false);
+        }
+      }
+      for (auto f : forest) {
+        for (auto sf : f.second) {
+          for (auto leaf : sf.getNdTrips()) {
+            // only one reference node
+            trips.push_back(leaf.second.front());
+          }
+        }
+      }
+    } else {
+      for (auto t : evalFeed->getTrips()) {
+        trips.push_back(t.second);
+      }
+    }
+
     LOG(DEBUG) << "Evaluating " << path << "...";
     size_t i = 0;
-    for (const auto& oldTrip : evalFeed->getTrips()) {
-      LOG(DEBUG) << "@ " << ++i << "/" << evalFeed->getTrips().size();
-      if (!mots->count(oldTrip.second->getRoute()->getType())) continue;
-      auto newTrip = feed.getTrips().get(oldTrip.first);
+    for (const auto& oldTrip : trips) {
+      LOG(DEBUG) << "@ " << ++i << "/" << trips.size();
+      if (!mots->count(oldTrip->getRoute()->getType())) continue;
+      auto newTrip = feed.getTrips().get(oldTrip->getId());
       if (!newTrip) {
-        LOG(ERROR) << "Trip #" << oldTrip.first << " not present in " << path
+        LOG(ERROR) << "Trip #" << oldTrip->getId() << " not present in " << path
                    << ", skipping...";
         continue;
       }
-      (*colls)[myFeed].add(oldTrip.second, oldTrip.second->getShape(), newTrip,
+      (*colls)[myFeed].add(oldTrip, oldTrip->getShape(), newTrip,
                            newTrip->getShape());
     }
   }
@@ -90,6 +134,7 @@ int main(int argc, char** argv) {
   bool summarize = false;
   bool json = false;
   bool avg = false;
+  bool unique = false;
 
   for (int i = 1; i < argc; i++) {
     std::string cur = argv[i];
@@ -106,6 +151,8 @@ int main(int argc, char** argv) {
       summarize = true;
     } else if (cur == "--json") {
       json = true;
+    } else if (cur == "--unique") {
+      unique = true;
     } else if (cur == "--avg") {
       avg = true;
     } else if (cur == "-f") {
@@ -169,8 +216,8 @@ int main(int argc, char** argv) {
 
   std::vector<std::thread> thrds(THREADS);
   for (auto& thr : thrds)
-    thr =
-        std::thread(&eval, &evlFeedPaths, &evalColls, &mots, &groundTruthFeed);
+    thr = std::thread(&eval, &evlFeedPaths, &evalColls, &mots, &groundTruthFeed,
+                      unique);
 
   for (auto& thr : thrds) thr.join();
 
@@ -188,9 +235,7 @@ int main(int argc, char** argv) {
     util::json::Dict jsonStats;
 
     if (evalColls.size() == 1) {
-      jsonStats = {
-          {"statistics", stats[evlFeedPaths[0]]
-      }};
+      jsonStats = {{"statistics", stats[evlFeedPaths[0]]}};
     } else {
       if (avg) {
         double count = evalColls.size();
@@ -206,13 +251,9 @@ int main(int argc, char** argv) {
           }
           avgStats[k] = sum / count;
         }
-        jsonStats = {
-            {"statistics", avgStats
-        }};
+        jsonStats = {{"statistics", avgStats}};
       } else {
-        jsonStats = {
-            {"statistics", stats
-        }};
+        jsonStats = {{"statistics", stats}};
       }
     }
 
diff --git a/src/util/geo/Geo.h b/src/util/geo/Geo.h
index 255b707..8b7f966 100644
--- a/src/util/geo/Geo.h
+++ b/src/util/geo/Geo.h
@@ -933,7 +933,7 @@ inline Line<T> lineFromWKT(std::string wkt) {
 template <typename T>
 inline std::string getWKT(const Point<T>& p) {
   std::stringstream ss;
-  ss << "POINT (" << p.getX() << " " << p.getY() << ")";
+  ss << std::fixed << std::setprecision(6) << "POINT (" << p.getX() << " " << p.getY() << ")";
   return ss.str();
 }
 
@@ -941,7 +941,7 @@ inline std::string getWKT(const Point<T>& p) {
 template <typename T>
 inline std::string getWKT(const std::vector<Point<T>>& p) {
   std::stringstream ss;
-  ss << "MULTIPOINT (";
+  ss << std::fixed << std::setprecision(6) << "MULTIPOINT (";
   for (size_t i = 0; i < p.size(); i++) {
     if (i) ss << ", ";
     ss << "(" << p[i].getX() << " " << p[i].getY() << ")";
@@ -954,7 +954,7 @@ inline std::string getWKT(const std::vector<Point<T>>& p) {
 template <typename T>
 inline std::string getWKT(const Line<T>& l) {
   std::stringstream ss;
-  ss << "LINESTRING (";
+  ss << std::fixed << std::setprecision(6) << "LINESTRING (";
   for (size_t i = 0; i < l.size(); i++) {
     if (i) ss << ", ";
     ss << l[i].getX() << " " << l[i].getY();
@@ -967,7 +967,7 @@ inline std::string getWKT(const Line<T>& l) {
 template <typename T>
 inline std::string getWKT(const std::vector<Line<T>>& ls) {
   std::stringstream ss;
-  ss << "MULTILINESTRING (";
+  ss << std::fixed << std::setprecision(6) << "MULTILINESTRING (";
 
   for (size_t j = 0; j < ls.size(); j++) {
     if (j) ss << ", ";