Skip to content

Commit 56181c0

Browse files
luca-sluca-s
committed
ENH: solver now reports more interesting statistics
1 parent 4aef108 commit 56181c0

File tree

7 files changed

+119
-105
lines changed

7 files changed

+119
-105
lines changed

‎apps/scrtdd/hdd/catalog.cpp

+7-7
Original file line numberDiff line numberDiff line change
@@ -306,12 +306,12 @@ Catalog::Catalog(const string& stationFile,
306306
{
307307
ev.relocInfo.isRelocated = true;
308308
ev.relocInfo.numNeighbours = std::stoul(row.at("numNeighbours"));
309+
ev.relocInfo.usedP = std::stoul(row.at("usedP"));
310+
ev.relocInfo.usedS = std::stoul(row.at("usedS"));
309311
ev.relocInfo.numCCp = std::stoul(row.at("numCCp"));
310312
ev.relocInfo.numCCs = std::stoul(row.at("numCCs"));
311313
ev.relocInfo.numCTp = std::stoul(row.at("numCTp"));
312314
ev.relocInfo.numCTs = std::stoul(row.at("numCTs"));
313-
ev.relocInfo.numObservs = std::stoul(row.at("numObservs"));
314-
ev.relocInfo.numFinalObservs = std::stoul(row.at("numFinalObservs"));
315315
ev.relocInfo.meanObsWeight = std::stod(row.at("meanObsWeight"));
316316
ev.relocInfo.meanFinalObsWeight = std::stod(row.at("meanFinalObsWeight"));
317317
}
@@ -342,7 +342,7 @@ Catalog::Catalog(const string& stationFile,
342342
ph.relocInfo.finalWeight = std::stod(row.at("finalWeight"));
343343
ph.relocInfo.residual = std::stod(row.at("residual"));
344344
ph.relocInfo.numObservs = std::stoul(row.at("numObservs"));
345-
ph.relocInfo.numFinalObservs = std::stoul(row.at("numFinalObservs"));
345+
ph.relocInfo.numXcorrObservs = std::stoul(row.at("numXcorrObservs"));
346346
ph.relocInfo.meanObsWeight = std::stod(row.at("meanObsWeight"));
347347
ph.relocInfo.meanFinalObsWeight = std::stod(row.at("meanFinalObsWeight"));
348348
}
@@ -754,7 +754,7 @@ void Catalog::writeToFile(string eventFile, string phaseFile, string stationFile
754754
stringstream evStreamReloc;
755755

756756
evStreamNoReloc << "id,isotime,latitude,longitude,depth,magnitude,rms";
757-
evStreamReloc << evStreamNoReloc.str() << ",relocated,numNeighbours,numCCp,numCCs,numCTp,numCTs,numObservs,numFinalObservs,meanObsWeight,meanFinalObsWeight" << endl;
757+
evStreamReloc << evStreamNoReloc.str() << ",relocated,numNeighbours,numPhaseP,numPhaseS,numCCp,numCCs,numCTp,numCTs,meanObsWeight,meanFinalObsWeight" << endl;
758758
evStreamNoReloc << endl;
759759

760760
bool relocInfo = false;
@@ -779,9 +779,9 @@ void Catalog::writeToFile(string eventFile, string phaseFile, string stationFile
779779
relocInfo = true;
780780
evStreamReloc << stringify(",true,%u,%u,%u,%u,%u,%u,%u,%.2f,%.2f",
781781
ev.relocInfo.numNeighbours,
782+
ev.relocInfo.usedP, ev.relocInfo.usedS,
782783
ev.relocInfo.numCCp, ev.relocInfo.numCCs,
783784
ev.relocInfo.numCTp, ev.relocInfo.numCTs,
784-
ev.relocInfo.numObservs, ev.relocInfo.numFinalObservs,
785785
ev.relocInfo.meanObsWeight, ev.relocInfo.meanFinalObsWeight);
786786
}
787787
evStreamReloc << endl;
@@ -798,7 +798,7 @@ void Catalog::writeToFile(string eventFile, string phaseFile, string stationFile
798798
phStream << "eventId,stationId,isotime,lowerUncertainty,upperUncertainty,type,networkCode,stationCode,locationCode,channelCode,evalMode";
799799
if (relocInfo)
800800
{
801-
phStream << ",usedInReloc,residual,initialWeight,finalWeight,numObservs,numFinalObservs,meanObsWeight,meanFinalObsWeight";
801+
phStream << ",usedInReloc,residual,initialWeight,finalWeight,numObservs,numXcorrObservs,meanObsWeight,meanFinalObsWeight";
802802
}
803803
phStream << endl;
804804

@@ -824,7 +824,7 @@ void Catalog::writeToFile(string eventFile, string phaseFile, string stationFile
824824
phStream << stringify(",true,%.3f,%.2f,%.2f,%u,%u,%.2f,%.2f",
825825
ph.relocInfo.residual,
826826
ph.procInfo.weight, ph.relocInfo.finalWeight,
827-
ph.relocInfo.numObservs, ph.relocInfo.numFinalObservs,
827+
ph.relocInfo.numObservs, ph.relocInfo.numXcorrObservs,
828828
ph.relocInfo.meanObsWeight, ph.relocInfo.meanFinalObsWeight);
829829
}
830830
}

‎apps/scrtdd/hdd/catalog.h

+3-3
Original file line numberDiff line numberDiff line change
@@ -80,12 +80,12 @@ class Catalog : public Core::BaseObject {
8080
struct {
8181
bool isRelocated = false;
8282
unsigned numNeighbours;
83+
unsigned usedP;
84+
unsigned usedS;
8385
unsigned numCCp;
8486
unsigned numCCs;
8587
unsigned numCTp;
8688
unsigned numCTs;
87-
unsigned numObservs;
88-
unsigned numFinalObservs;
8989
double meanObsWeight;
9090
double meanFinalObsWeight;
9191
} relocInfo;
@@ -138,7 +138,7 @@ class Catalog : public Core::BaseObject {
138138
double residual;
139139
double finalWeight;
140140
unsigned numObservs;
141-
unsigned numFinalObservs;
141+
unsigned numXcorrObservs;
142142
double meanObsWeight;
143143
double meanFinalObsWeight;
144144
} relocInfo;

‎apps/scrtdd/hdd/hypodd.cpp

+68-68
Original file line numberDiff line numberDiff line change
@@ -465,8 +465,7 @@ HypoDD::relocate(CatalogPtr& catalog,
465465
//
466466
for (const NeighboursPtr& neighbours : neighbourCats)
467467
{
468-
catalog = addObservations(solver, catalog, neighbours, keepNeighboursFixed,
469-
xcorr, obsparams);
468+
addObservations(solver, catalog, neighbours, keepNeighboursFixed, xcorr, obsparams);
470469
}
471470
obsparams.addToSolver(solver);
472471

@@ -493,6 +492,12 @@ HypoDD::relocate(CatalogPtr& catalog,
493492
//
494493
// solve the system
495494
//
495+
SEISCOMP_INFO("Solving iteration %lu num events %u (dampingFactor %.2f "
496+
"downWeightingByResidual %.2f meanShiftConstrainWeight %.2f,%.2f,%.2f,%.2f)",
497+
iteration, neighbourCats.size(), dampingFactor, downWeightingByResidual,
498+
meanLonShiftConstraint, meanLatShiftConstraint,
499+
meanDepthShiftConstraint, meanTTShiftConstraint);
500+
496501
try {
497502
solver.solve(_cfg.solver.solverIterations, dampingFactor,
498503
downWeightingByResidual, meanLonShiftConstraint,
@@ -508,22 +513,6 @@ HypoDD::relocate(CatalogPtr& catalog,
508513

509514
// update event parameters
510515
catalog = updateRelocatedEvents(solver, catalog, neighbourCats, obsparams);
511-
512-
vector<double> rms;
513-
for (const NeighboursPtr& neighbours : neighbourCats)
514-
{
515-
const Event& ev = catalog->getEvents().at(neighbours->refEvId);
516-
if ( ev.relocInfo.isRelocated ) rms.push_back(ev.rms);
517-
}
518-
SEISCOMP_INFO("Iteration %u dampingFactor %.2f downWeightingByResidual %.2f "
519-
"meanShiftConstrainWeight %.2f,%.2f,%.2f,%.2f",
520-
iteration, dampingFactor, downWeightingByResidual,
521-
meanLonShiftConstraint, meanLatShiftConstraint,
522-
meanDepthShiftConstraint, meanTTShiftConstraint);
523-
const double median = computeMedian(rms);
524-
const double MAD = computeMedianAbsoluteDeviation(rms, median);
525-
SEISCOMP_INFO("Events Rms: median %.4f median absolute deviation %.4f",
526-
median*1000, MAD*1000);
527516
}
528517

529518
// build the relocated catalog from the results of relocations
@@ -547,13 +536,15 @@ string HypoDD::relocationReport(const CatalogCPtr& relocatedEv)
547536
if ( ! event.relocInfo.isRelocated )
548537
return "Event not relocated";
549538

550-
return stringify("Rms residual %.3f [sec]. Neighboring events %u "
551-
"Num observations: initial %u (avg. weight %.2f) final %u (avg. weight %.2f). "
552-
"Num cross-correlated observations: P phases %u S phases %u. "
539+
return stringify("Rms residual %.3f [sec]. Neighboring events %u used P %u used S %u."
540+
"Num observations: %u (avg. a priory weight %.2f final weight %.2f). "
541+
"Num xcross observations: P phases %u S phases %u. "
553542
"Num absolute TT diff observations: P phases %u S phases %u",
554-
event.rms, event.relocInfo.numNeighbours,
555-
event.relocInfo.numObservs, event.relocInfo.meanObsWeight,
556-
event.relocInfo.numFinalObservs, event.relocInfo.meanFinalObsWeight,
543+
event.rms, event.relocInfo.numNeighbours,
544+
event.relocInfo.usedP, event.relocInfo.usedS,
545+
(event.relocInfo.numCCp+event.relocInfo.numCCs+
546+
event.relocInfo.numCTp+event.relocInfo.numCTs),
547+
event.relocInfo.meanObsWeight, event.relocInfo.meanFinalObsWeight,
557548
event.relocInfo.numCCp, event.relocInfo.numCCs,
558549
event.relocInfo.numCTp, event.relocInfo.numCTs);
559550
}
@@ -564,17 +555,13 @@ string HypoDD::relocationReport(const CatalogCPtr& relocatedEv)
564555
* differential travel times from cross correlation for pairs of
565556
* earthquakes
566557
*/
567-
CatalogPtr
558+
void
568559
HypoDD::addObservations(Solver& solver, const CatalogCPtr& catalog, const NeighboursPtr& neighbours,
569560
bool keepNeighboursFixed, const XCorrCache& xcorr,
570561
ObservationParams& obsparams ) const
571562
{
572563
// copy event because we'll update it
573-
Event refEv = catalog->getEvents().at(neighbours->refEvId);
574-
refEv.relocInfo.numCCp = 0;
575-
refEv.relocInfo.numCCs = 0;
576-
refEv.relocInfo.numCTp = 0;
577-
refEv.relocInfo.numCTs = 0;
564+
const Event& refEv = catalog->getEvents().at(neighbours->refEvId);
578565

579566
//
580567
// loop through reference event phases
@@ -634,6 +621,7 @@ HypoDD::addObservations(Solver& solver, const CatalogCPtr& catalog, const Neighb
634621
double weight = _cfg.solver.usePickUncertainty
635622
? (refPhase.procInfo.weight + phase.procInfo.weight) / 2.0
636623
: 1.0;
624+
bool isXcorr = false;
637625
//
638626
// Check if we have xcorr results for current event/refEvent pair at station/phase
639627
// and use those instead
@@ -644,25 +632,17 @@ HypoDD::addObservations(Solver& solver, const CatalogCPtr& catalog, const Neighb
644632
const auto& xcdata = xcorr.get(refEv.id, event.id, refPhase.stationId, refPhase.procInfo.type);
645633
diffTime -= xcdata.lag;
646634
weight *= _cfg.solver.xcorrObsWeight;
647-
if (refPhase.procInfo.type == Phase::Type::P) refEv.relocInfo.numCCp++;
648-
if (refPhase.procInfo.type == Phase::Type::S) refEv.relocInfo.numCCs++;
635+
isXcorr = true;
649636
}
650637
else
651638
{
652639
weight *= _cfg.solver.absTTDiffObsWeight;
653-
if (refPhase.procInfo.type == Phase::Type::P) refEv.relocInfo.numCTp++;
654-
if (refPhase.procInfo.type == Phase::Type::S) refEv.relocInfo.numCTs++;
655640
}
656641

657642
solver.addObservation(refEv.id, event.id, refPhase.stationId, phaseTypeAsChar,
658-
diffTime, weight, true, !keepNeighboursFixed);
643+
diffTime, weight, true, !keepNeighboursFixed, isXcorr);
659644
}
660645
}
661-
662-
// save back the computed refEv.relocInfo.numCC/CT P/S information
663-
CatalogPtr returnCat(new Catalog(*catalog) );
664-
returnCat->updateEvent(refEv);
665-
return returnCat;
666646
}
667647

668648
void
@@ -710,6 +690,7 @@ HypoDD::updateRelocatedEvents(const Solver& solver,
710690
unordered_map<string,Station> stations = catalog->getStations();
711691
map<unsigned,Event> events = catalog->getEvents();
712692
unordered_multimap<unsigned,Phase> phases = catalog->getPhases();
693+
unsigned relocatedEvs = 0;
713694

714695
for (const NeighboursPtr& neighbours : neighbourCats)
715696
{
@@ -728,19 +709,26 @@ HypoDD::updateRelocatedEvents(const Solver& solver,
728709
continue;
729710
}
730711

712+
relocatedEvs++;
713+
731714
event.latitude += deltaLat;
732715
event.longitude += deltaLon;
733716
event.depth += deltaDepth;
734717
event.time += Core::TimeSpan(deltaTT);
735718
event.relocInfo.isRelocated = true;
736719
event.relocInfo.numNeighbours = neighbours->numNeighbours;
737720

738-
event.rms = 0.;
739-
event.relocInfo.numObservs = 0;
740-
event.relocInfo.numFinalObservs = 0;
721+
event.relocInfo.usedP = 0;
722+
event.relocInfo.usedS = 0;
723+
event.relocInfo.numCCp = 0;
724+
event.relocInfo.numCCs = 0;
725+
event.relocInfo.numCTp = 0;
726+
event.relocInfo.numCTs = 0;
741727
event.relocInfo.meanObsWeight = 0;
742728
event.relocInfo.meanFinalObsWeight = 0;
729+
event.rms = 0.;
743730

731+
unsigned rmsCount = 0;
744732
unsigned pCount = 0;
745733
auto eqlrng = phases.equal_range(event.id);
746734
for (auto it = eqlrng.first; it != eqlrng.second; ++it)
@@ -751,52 +739,64 @@ HypoDD::updateRelocatedEvents(const Solver& solver,
751739

752740
phase.relocInfo.isRelocated = false;
753741

754-
unsigned startingObservations, finalObservations;
755-
double totalAPrioriWeight, totalFinalWeight;
742+
unsigned startingObservations, startingXcorrObservations, totalFinalObservations;
743+
double meanAPrioriWeight, meanFinalWeight;
756744
if ( ! solver.getObservationParamsChanges(event.id, station.id, phaseTypeAsChar,
757-
startingObservations, finalObservations,
758-
totalAPrioriWeight, totalFinalWeight) )
745+
startingObservations,
746+
startingXcorrObservations,
747+
totalFinalObservations,
748+
meanAPrioriWeight,
749+
meanFinalWeight) )
759750
{
760751
continue;
761752
}
762753

763-
if ( finalObservations == 0 )
764-
continue;
754+
event.relocInfo.meanObsWeight += meanAPrioriWeight;
755+
event.relocInfo.meanFinalObsWeight += meanFinalWeight;
756+
++pCount;
757+
758+
if ( meanFinalWeight == 0 || totalFinalObservations == 0 )
759+
continue;
765760

766761
phase.relocInfo.isRelocated = true;
767-
phase.relocInfo.finalWeight = phase.procInfo.weight * totalFinalWeight / totalAPrioriWeight;
768-
phase.relocInfo.numObservs = startingObservations;
769-
phase.relocInfo.numFinalObservs = finalObservations;
770-
phase.relocInfo.meanObsWeight = totalAPrioriWeight / startingObservations;
771-
phase.relocInfo.meanFinalObsWeight = totalFinalWeight / finalObservations;
762+
phase.relocInfo.finalWeight = phase.procInfo.weight * meanFinalWeight / meanAPrioriWeight;
763+
phase.relocInfo.numObservs = startingObservations;
764+
phase.relocInfo.numXcorrObservs = startingXcorrObservations;
765+
phase.relocInfo.meanObsWeight = meanAPrioriWeight;
766+
phase.relocInfo.meanFinalObsWeight = meanFinalWeight;
772767

773768
try {
774769
obsparams.add(_ttt, event, station, phaseTypeAsChar);
775-
} catch ( exception &e ) {
776-
SEISCOMP_DEBUG("Skipping observation parameter (ev %u sta %s phase %c): %s",
777-
event.id, station.id.c_str(), phaseTypeAsChar, e.what());
778-
continue;
779-
}
780-
781-
double travelTime = obsparams.get(event.id, station.id, phaseTypeAsChar).travelTime;
782-
phase.relocInfo.residual = travelTime - (phase.time - event.time).length();
770+
double travelTime = obsparams.get(event.id, station.id, phaseTypeAsChar).travelTime;
771+
phase.relocInfo.residual = travelTime - (phase.time - event.time).length();
772+
rmsCount++;
773+
} catch ( exception &e ) { }
783774

784775
event.rms += (phase.relocInfo.residual * phase.relocInfo.residual);
785-
event.relocInfo.numObservs += phase.relocInfo.numObservs;
786-
event.relocInfo.numFinalObservs += phase.relocInfo.numFinalObservs;
787-
event.relocInfo.meanObsWeight += phase.relocInfo.meanObsWeight;
788-
event.relocInfo.meanFinalObsWeight += phase.relocInfo.meanFinalObsWeight;
789-
++pCount;
776+
if ( phase.procInfo.type == Phase::Type::P )
777+
{
778+
event.relocInfo.usedP++;
779+
event.relocInfo.numCCp += phase.relocInfo.numXcorrObservs;
780+
event.relocInfo.numCTp += phase.relocInfo.numObservs;
781+
}
782+
if ( phase.procInfo.type == Phase::Type::S )
783+
{
784+
event.relocInfo.usedS++;
785+
event.relocInfo.numCCs += phase.relocInfo.numXcorrObservs;
786+
event.relocInfo.numCTs += phase.relocInfo.numObservs;
787+
}
790788
}
791789

790+
if ( rmsCount > 0 ) event.rms = std::sqrt(event.rms / rmsCount);
792791
if ( pCount > 0 )
793792
{
794-
event.rms = std::sqrt(event.rms / pCount);
795793
event.relocInfo.meanObsWeight /= pCount;
796794
event.relocInfo.meanFinalObsWeight /= pCount;
797795
}
798796
}
799797

798+
SEISCOMP_INFO("Successfully relocated %u events", relocatedEvs);
799+
800800
return new Catalog(stations, events, phases);
801801
}
802802

‎apps/scrtdd/hdd/hypodd.h

+4-4
Original file line numberDiff line numberDiff line change
@@ -192,10 +192,10 @@ class HypoDD : public Core::BaseObject {
192192
std::unordered_map<std::string,Entry> _entries;
193193
};
194194

195-
CatalogPtr addObservations(Solver& solver, const CatalogCPtr& catalog,
196-
const NeighboursPtr& neighbours,
197-
bool keepNeighboursFixed, const XCorrCache& xcorr,
198-
ObservationParams& obsparams ) const;
195+
void addObservations(Solver& solver, const CatalogCPtr& catalog,
196+
const NeighboursPtr& neighbours,
197+
bool keepNeighboursFixed, const XCorrCache& xcorr,
198+
ObservationParams& obsparams ) const;
199199

200200
CatalogPtr updateRelocatedEvents(const Solver& solver,
201201
const CatalogCPtr& catalog,

0 commit comments

Comments
 (0)