ENH: multi-event relocation is now performed per cluster (vs all events in one)

Instead of relocating all events in one single double-difference system
each cluster is relocated independently. This results in improved
solutions since the parameters of the solver are cluster dependent
(e.g. downweighting by residuals)
Also each cluster covers a smaller area, which helps in avoiding
breaking the solver assumption of an euclidean geometry (the lat/lon/depth
coordinates are converted internally to an euclidean system by the solver)

Author: luca-s

apps/scrtdd/hdd/clustering.cpp

@@ -346,8 +346,7 @@ selectNeighbouringEvents(const CatalogCPtr& catalog,
 }
 
 
-
-std::list<NeighboursPtr>
+deque< list<NeighboursPtr> >
 selectNeighbouringEventsCatalog(const CatalogCPtr& catalog,
                                 double minPhaseWeight,
                                 double minESdist,
@@ -364,7 +363,7 @@ selectNeighbouringEventsCatalog(const CatalogCPtr& catalog,
     SEISCOMP_INFO("Selecting Catalog Neighbouring Events ");
 
     // neighbours for each event
-    list<NeighboursPtr> neighboursByEvent;
+    list<NeighboursPtr> neighboursList;
 
     // for each event find the neighbours
     CatalogPtr validCatalog = new Catalog(*catalog);
@@ -408,7 +407,7 @@ selectNeighbouringEventsCatalog(const CatalogCPtr& catalog,
         // add newly computed neighbors catalogs to previous ones
         for ( NeighboursPtr& neighbours : newNeighbourCats )
         {
-            neighboursByEvent.push_back( neighbours );
+            neighboursList.push_back( neighbours );
             // make sure we won't recompute what has been already done
             todoEvents.remove( neighbours->refEvId );
         }
@@ -420,7 +419,7 @@ selectNeighbouringEventsCatalog(const CatalogCPtr& catalog,
             redo = false;
             list<NeighboursPtr> validNeighbourCats;
 
-            for ( NeighboursPtr& neighbours : neighboursByEvent )
+            for ( NeighboursPtr& neighbours : neighboursList )
             {
                 bool currCatInvalid = false;
                 for (unsigned removedEventId : removedEvents)
@@ -442,38 +441,107 @@ selectNeighbouringEventsCatalog(const CatalogCPtr& catalog,
                 validNeighbourCats.push_back( neighbours );
             }
 
-            neighboursByEvent.clear();
-            neighboursByEvent = validNeighbourCats;
+            neighboursList.clear();
+            neighboursList = validNeighbourCats;
 
         } while( redo );
     }
 
-    // We don't want to report the same pairs multiple times
-    // when creating double-difference observations. So we'll
-    // remove the pairs that appeared in previous catalogs from
-    // the following catalogs
-    std::unordered_multimap<unsigned,unsigned> existingPairs;
+    return clusterizeNeighbouringEvents(neighboursList);
+}
+
+
+/*
+ * Organize the neighbours by not connected clusters
+ * Also, we don't want to report the same pair multiple times
+ * (e.g. ev1-ev2 and ev2-ev1) since we only want one observation
+ * for pair when creating the double-difference observations
+ * system
+ */
+deque< list<NeighboursPtr> >
+clusterizeNeighbouringEvents(const list<NeighboursPtr>& neighboursList)
+{
+    map<unsigned, list<NeighboursPtr> > clusters;
+
+    unordered_map<unsigned, unsigned> clusterIdByEvent; // event id, cluster id
+
+    unordered_map<unsigned, NeighboursPtr> neighboursByEvent; // key event id
+    for ( const NeighboursPtr& neighbours : neighboursList )
+        neighboursByEvent[ neighbours->refEvId ] = neighbours;
 
-    for ( NeighboursPtr& neighbours : neighboursByEvent )
+    while ( ! neighboursByEvent.empty() )
     {
-        unsigned currEventId = neighbours->refEvId;
+        // keep track of event pairs found (for dropping identical pairs)
+        unordered_multimap<unsigned,unsigned> discoveredPairs;
 
-        // remove from currrent catalog the existing pairs
-        auto eqlrng = existingPairs.equal_range(currEventId );
-        for (auto existingPair = eqlrng.first; existingPair != eqlrng.second; existingPair++)
+        // start traversal with first unseen event neighbours
+        unordered_set<unsigned> clusterEvs( { neighboursByEvent.begin()->first } );
+
+        list<NeighboursPtr> currentCluster;
+        unordered_set<unsigned> connectedClusters;
+
+        while ( ! clusterEvs.empty() ) // when empty the cluster is fully built
         {
-            neighbours->ids.erase( existingPair->second );
-            neighbours->phases.erase( existingPair->second );
+            unsigned currentEv = *clusterEvs.begin();
+            clusterEvs.erase(currentEv);
+
+            // keep track of clusters connected to the current one
+            if ( clusterIdByEvent.find(currentEv) != clusterIdByEvent.end() )
+                connectedClusters.insert( clusterIdByEvent.at(currentEv) );
+
+            const auto& neighboursByEventIt = neighboursByEvent.find( currentEv );
+
+            // skip already processed events
+            if ( neighboursByEventIt == neighboursByEvent.end() )
+                continue;
+
+            NeighboursPtr neighbours = neighboursByEventIt->second;
+            neighboursByEvent.erase( neighbours->refEvId );
+
+            // update the set for the traversal of this cluster
+            for ( unsigned neighEvId : neighbours->ids )
+                clusterEvs.insert( neighEvId );
+
+            // remove from current neighbours the pairs that appeared previously
+            // in this cluster
+            auto eqlrng = discoveredPairs.equal_range( neighbours->refEvId );
+            for (auto existingPair = eqlrng.first; existingPair != eqlrng.second; existingPair++)
+            {
+                neighbours->ids.erase( existingPair->second );
+                neighbours->phases.erase( existingPair->second );
+            }
+
+            // keep track of new event pairs for following iterations and
+            // update the queue for the breadth-first traversal
+            for ( unsigned neighEvId : neighbours->ids )
+                discoveredPairs.emplace(neighEvId, neighbours->refEvId);
+
+            // populate current cluster
+            currentCluster.push_back( neighbours );
         }
 
-        // remove current pairs from following catalogs
-        for ( unsigned neighEvId : neighbours->ids )
+        if ( ! connectedClusters.empty() )
         {
-            existingPairs.emplace(neighEvId, currEventId);
+            // merge all connected clusters to the current one
+            for ( unsigned clusterId : connectedClusters )
+            {
+                currentCluster.splice(currentCluster.end(), clusters[clusterId] );
+                clusters.erase(clusterId);
+            }
         }
+
+        unsigned maxKey = clusters.empty() ? 0 : clusters.rbegin()->first;
+        unsigned newClusterId = maxKey + 1;
+
+        clusters[newClusterId] = currentCluster;
+
+        for ( const NeighboursPtr& n : currentCluster )
+            clusterIdByEvent[ n->refEvId ] = newClusterId;
     }
 
-    return neighboursByEvent;
+    deque< list<NeighboursPtr> > returnClusters;
+    for (auto& kv : clusters ) returnClusters.push_back(kv.second);
+    return returnClusters;
 }
 
 

apps/scrtdd/hdd/clustering.h

@@ -23,6 +23,8 @@
 #include <unordered_map>
 #include <unordered_set>
 #include <set>
+#include <list>
+#include <deque>
 
 
 namespace Seiscomp {
@@ -104,7 +106,7 @@ selectNeighbouringEvents(const CatalogCPtr& catalog,
                          double maxEllipsoidSize=10,
                          bool keepUnmatched=false);
 
-std::list<NeighboursPtr>
+std::deque< std::list<NeighboursPtr> >
 selectNeighbouringEventsCatalog(const CatalogCPtr& catalog,
                                 double minPhaseWeight,
                                 double minESdis,
@@ -118,6 +120,9 @@ selectNeighbouringEventsCatalog(const CatalogCPtr& catalog,
                                 double maxEllipsoidSize,
                                 bool keepUnmatched);
 
+std::deque< std::list<NeighboursPtr> >
+clusterizeNeighbouringEvents(const std::list<NeighboursPtr>& neighboursList);
+
 }
 }
 

apps/scrtdd/hdd/hypodd.cpp

@@ -228,7 +228,7 @@ CatalogPtr HypoDD::relocateCatalog()
     }
 
     // Find Neighbouring Events in the catalog
-    list<NeighboursPtr> neighbourCats = selectNeighbouringEventsCatalog(
+    deque< list<NeighboursPtr> > clusters = selectNeighbouringEventsCatalog(
         catToReloc, _cfg.ddObservations2.minWeight,
         _cfg.ddObservations2.minESdist, _cfg.ddObservations2.maxESdist,
         _cfg.ddObservations2.minEStoIEratio, _cfg.ddObservations2.minDTperEvt,
@@ -237,28 +237,58 @@ CatalogPtr HypoDD::relocateCatalog()
         _cfg.ddObservations2.maxEllipsoidSize, true
     );
 
+    SEISCOMP_INFO("Found %lu event clusters", clusters.size() );
+
     // write catalog for debugging purpose
     if ( ! _workingDirCleanup )
     {
         catToReloc->writeToFile(
             (boost::filesystem::path(catalogWorkingDir)/"starting-event.csv").string(),
             (boost::filesystem::path(catalogWorkingDir)/"starting-phase.csv").string(),
             (boost::filesystem::path(catalogWorkingDir)/"starting-station.csv").string() );
-        CatalogPtr catToDump( new Catalog() );
-        for (const auto& n : neighbourCats)
-            catToDump->add(n->refEvId, *catToReloc, true);
-        catToDump->writeToFile(
-            (boost::filesystem::path(catalogWorkingDir)/"selected-event.csv").string(),
-            (boost::filesystem::path(catalogWorkingDir)/"selected-phase.csv").string(),
-            (boost::filesystem::path(catalogWorkingDir)/"selected-station.csv").string() ); 
-    } 
+    }
+
+    //
+    // Relocate one cluster a time
+    //
+    CatalogPtr relocatedCatalog( new Catalog() );
+
+    for (unsigned clusterId = 0; clusterId < clusters.size(); clusterId++ )
+    {
+        const list<NeighboursPtr>& neighCluster = clusters[clusterId];
+
+        SEISCOMP_INFO("Relocating cluster %u (%lu events)", clusterId+1, neighCluster.size() );
+
+        if ( ! _workingDirCleanup )
+        {
+            CatalogPtr catToDump( new Catalog() );
+            for (const NeighboursPtr& n : neighCluster )
+                catToDump->add(n->refEvId, *catToReloc, true);
+            string prefix = "cluster-" + to_string(clusterId + 1);
+            catToDump->writeToFile(
+                (boost::filesystem::path(catalogWorkingDir)/(prefix+"-event.csv")).string(),
+                (boost::filesystem::path(catalogWorkingDir)/(prefix+"-phase.csv")).string(),
+                (boost::filesystem::path(catalogWorkingDir)/(prefix+"-station.csv")).string() ); 
+        }
 
-    // Perform cross correlation, which also detects picks around theoretical
-    // arrival times. The catalog will be updated with those theoretical phases 
-    const XCorrCache xcorr = buildXCorrCache(catToReloc, neighbourCats, _useArtificialPhases);
+        // Perform cross correlation, which also detects picks around theoretical
+        // arrival times. The catalog will be updated with those theoretical phases 
+        const XCorrCache xcorr = buildXCorrCache(catToReloc, neighCluster, _useArtificialPhases);
 
-    // The actual relocation
-    CatalogPtr relocatedCatalog = relocate(catToReloc, neighbourCats, false, xcorr);
+        // The actual relocation
+        CatalogPtr relocatedCluster = relocate(catToReloc, neighCluster, false, xcorr);
+
+        relocatedCatalog->add(*relocatedCluster, true);
+
+        if ( ! _workingDirCleanup )
+        {
+            string prefix = "relocated-cluster-" + to_string(clusterId + 1);
+            relocatedCluster->writeToFile(
+                (boost::filesystem::path(catalogWorkingDir)/(prefix+"-event.csv")).string(),
+                (boost::filesystem::path(catalogWorkingDir)/(prefix+"-phase.csv")).string(),
+                (boost::filesystem::path(catalogWorkingDir)/(prefix+"-station.csv")).string() ); 
+        }
+    }
 
     // write catalog for debugging purpose
     if ( ! _workingDirCleanup )
@@ -448,12 +478,10 @@ HypoDD::relocateEventSingleStep(const CatalogCPtr& evToRelocateCat,
 
 CatalogPtr
 HypoDD::relocate(CatalogPtr& catalog,
-                 const std::list<NeighboursPtr>& neighbourCats,
+                 const std::list<NeighboursPtr>& neighCluster,
                  bool keepNeighboursFixed,
                  const XCorrCache& xcorr) const
 {
-    // Create a solver and then add observations
-    Solver solver(_cfg.solver.type);
     ObservationParams obsparams;
 
     for ( unsigned iteration=0; iteration < _cfg.solver.algoIterations; iteration++ )
@@ -485,18 +513,18 @@ HypoDD::relocate(CatalogPtr& catalog,
         SEISCOMP_INFO("Solving iteration %u num events %lu (observ. weight TTdiff %.2f xcorr %.2f "
                       "dampingFactor %.2f downWeightingByResidual %.2f "
                       "meanShiftConstrainWeight %.2f,%.2f,%.2f,%.2f)",
-                      iteration, neighbourCats.size(), absTTDiffObsWeight, xcorrObsWeight,
+                      iteration, neighCluster.size(), absTTDiffObsWeight, xcorrObsWeight,
                       dampingFactor, downWeightingByResidual,
                       meanLonShiftConstraint, meanLatShiftConstraint,
                       meanDepthShiftConstraint, meanTTShiftConstraint); 
 
-
-        solver.reset();
+        // Create a solver and then add observations
+        Solver solver(_cfg.solver.type);
 
         //
         // Add absolute travel time/xcorr differences to the solver (the observations)
         //
-        for (const NeighboursPtr& neighbours : neighbourCats)
+        for (const NeighboursPtr& neighbours : neighCluster)
         {
             addObservations(solver, absTTDiffObsWeight, xcorrObsWeight, catalog,
                             neighbours, keepNeighboursFixed, xcorr, obsparams);
@@ -521,12 +549,12 @@ HypoDD::relocate(CatalogPtr& catalog,
         obsparams = ObservationParams();
 
         // update event parameters
-        catalog = updateRelocatedEvents(solver, catalog, neighbourCats, obsparams);
+        catalog = updateRelocatedEvents(solver, catalog, neighCluster, obsparams);
     }
 
     // build the relocated catalog from the results of relocations
     CatalogPtr relocatedCatalog( new Catalog() );
-    for (const NeighboursPtr& neighbours : neighbourCats)
+    for (const NeighboursPtr& neighbours : neighCluster)
     {
         const Event& event = catalog->getEvents().at(neighbours->refEvId);
         if ( event.relocInfo.isRelocated )
@@ -695,15 +723,15 @@ HypoDD::ObservationParams::addToSolver(Solver& solver) const
 CatalogPtr
 HypoDD::updateRelocatedEvents(const Solver& solver,
                               const CatalogCPtr& catalog,
-                              const std::list<NeighboursPtr>& neighbourCats,
+                              const std::list<NeighboursPtr>& neighCluster,
                               ObservationParams& obsparams ) const
 {
     unordered_map<string,Station> stations = catalog->getStations();
     map<unsigned,Event> events = catalog->getEvents();
     unordered_multimap<unsigned,Phase> phases = catalog->getPhases();
     unsigned relocatedEvs = 0;
 
-    for (const NeighboursPtr& neighbours : neighbourCats)
+    for (const NeighboursPtr& neighbours : neighCluster)
     {
         Event& event = events.at(neighbours->refEvId);
         event.relocInfo.isRelocated = false;
@@ -1014,14 +1042,14 @@ HypoDD::createThoreticalPhase(const Station& station,
 
 XCorrCache
 HypoDD::buildXCorrCache(CatalogPtr& catalog,
-                        const std::list<NeighboursPtr>& neighbourCats,
+                        const std::list<NeighboursPtr>& neighCluster,
                         bool computeTheoreticalPhases)
 {
     XCorrCache xcorr;
     _counters = {0};
     _wf->resetCounters();
 
-    for (const NeighboursPtr& neighbours : neighbourCats)
+    for (const NeighboursPtr& neighbours : neighCluster)
     {
         const Event& refEv = catalog->getEvents().at(neighbours->refEvId);