Merge pull request #70 from jpatton-USGS/publish-fix

Publish fix

cmake/cppcheck.cmake

@@ -31,6 +31,7 @@ if(RUN_CPPCHECK)
       --verbose
       --suppress=nullPointerRedundantCheck
       --suppress=constStatement
+      --suppress=compareBoolExpressionWithInt
       --error-exitcode=1
       ${SRCS} ${HDRS}
       COMMENT "Running cppcheck" VERBATIM

examples/global_example/initialize.d

@@ -15,6 +15,8 @@
       "AssociationStandardDeviationCutoff": 5.0,
       "PruningStandardDeviationCutoff": 5.0,
       "NonLocatingPhaseCutoffFactor": 3.0,
+      "TeleseismicDistanceLimit": 30.0,
+      "TeleseismicPhaseCountThreshold": 30,
       "PickAffinityExponentialFactor": 2.5,
       "DistanceCutoffFactor": 5.0,
       "DistanceCutoffRatio": 0.8,
@@ -125,6 +127,13 @@
           "TravFile": "./global_example/P'P'df.trv",
           "UseForLocation": false,
           "PublishPhase": false
+      },
+      {
+          "PhaseName": "PKiKP",
+          "Assoc": [ 55, 118 ],
+          "TravFile": "./global_example/PKiKP.trv",
+          "UseForLocation": false,
+          "PublishPhase": false
       }
   ]
 }

examples/global_example/output.d

@@ -10,6 +10,10 @@
 	# Whether to publish an event on expiration from glass
 	"PublishOnExpiration":true,
 
+	# Whether to immediately publish a quake that hits 
+	# this bayes threshold.
+	"ImmediatePublicationThreshold": 30.0,
+
 	# the directory to write output to
 	"OutputDirectory":"./global_example/output",
 

glasscore/glasslib/include/HypoList.h

@@ -310,6 +310,13 @@ class CHypoList : public glass3::util::ThreadBaseClass {
 	 */
 	glass3::util::WorkState work() override;
 
+	// constants
+	/**
+	 * \brief The duration in seconds to search into the past for hypos matching
+	 * a pick
+	 */
+	static constexpr double k_nHypoSearchPastDurationForPick = 3600;
+
  private:
 	/**
 	 * \brief A HypoList function that updates the position of the given hypo
@@ -370,12 +377,6 @@ class CHypoList : public glass3::util::ThreadBaseClass {
 	mutable std::recursive_mutex m_HypoListMutex;
 
 	// constants
-	/**
-	 * \brief The duration in seconds to search into the past for hypos matching
-	 * a pick
-	 */
-	static const int k_nHypoSearchPastDurationForPick = 2400;
-
 	/**
 	 * \brief default maximum number of hypos allowed in CHypoList.
 	 */

glasscore/glasslib/include/PickList.h

@@ -23,6 +23,7 @@
 
 #include "Glass.h"
 #include "Pick.h"
+#include "HypoList.h"
 
 namespace glasscore {
 
@@ -152,11 +153,12 @@ class CPickList : public glass3::util::ThreadBaseClass {
 	 * \param hyp - A shared_ptr to a CHypo object containing the hypocenter
 	 * to attempt to associate to.
 	 * \param tWindow - A double value containing the window to search picks
-	 * from origin time in seconds, defaults to 2400.0
+	 * from origin time in seconds, defaults to 3600.0
 	 * \return Returns true if any picks were associated to the hypocenter,
 	 * false otherwise.
 	 */
-	bool scavenge(std::shared_ptr<CHypo> hyp, double tWindow = 2400.0);
+	bool scavenge(std::shared_ptr<CHypo> hyp,
+		double tWindow = CHypoList::k_nHypoSearchPastDurationForPick);
 
 	/**
 	 * \brief Get the CSiteList pointer used by this pick list for site lookups

glasscore/glasslib/src/Hypo.cpp

@@ -1067,9 +1067,9 @@ bool CHypo::canAssociate(std::shared_ptr<CPick> pick, double sigma,
 		}
 	}
 
-	// get residula, get useForLocations
-        bool useForLocations = true;
-        double tRes = calculateResidual(pick, &useForLocations);
+	// get residual, get useForLocations
+	bool useForLocations = true;
+	double tRes = calculateResidual(pick, &useForLocations);
 
         // give up if there's no valid residual
 	if (std::isnan(tRes) == true) {
@@ -1107,7 +1107,7 @@ bool CHypo::canAssociate(std::shared_ptr<CPick> pick, double sigma,
 			// for locations, this will cause issues some sort of way to
 			// keep track of whether a phase is teleseismic will need to
 			// be added to traveltime
-			return (false);
+			// return (false);
 		}
 	}
 
@@ -1587,6 +1587,8 @@ double CHypo::calculateResidual(std::shared_ptr<CPick> pick,
 	// init expected travel time
 	double tCal = traveltime::CTravelTime::k_dTravelTimeInvalid;
 
+	std::string phaseName = "??";
+
 	// check pick classification
 	// are we configured to check pick phase classification
 	if (CGlass::getPickPhaseClassificationThreshold() > 0) {
@@ -1602,17 +1604,20 @@ double CHypo::calculateResidual(std::shared_ptr<CPick> pick,
 			// the classified pick phase
 			tCal = m_pTravelTimeTables->T(&site->getGeo(),
 											pick->getClassifiedPhase());
+			phaseName = pick->getClassifiedPhase();
 		} else {
 			// no valid phase classification,
 			// compute expected travel time based on the pick site location and
 			// the observed travel time
 			tCal = m_pTravelTimeTables->T(&site->getGeo(), tObs);
+			phaseName = m_pTravelTimeTables->m_sPhase;
 		}
 	} else {
 		// not checking phase classification
 		// compute expected travel time based on the pick site location and
 		// the observed travel time
 		tCal = m_pTravelTimeTables->T(&site->getGeo(), tObs);
+		phaseName = m_pTravelTimeTables->m_sPhase;
 	}
 
 	// Check if pick has an invalid travel time,
@@ -1631,6 +1636,24 @@ double CHypo::calculateResidual(std::shared_ptr<CPick> pick,
 
 	// compute residual from observed and calculated travel times
 	double tRes = tObs - tCal;
+
+	// set up geo for distance calculations
+	/* glass3::util::Geo geo;
+	geo.setGeographic(m_dLatitude, m_dLongitude,
+		glass3::util::Geo::k_EarthRadiusKm - m_dDepth);
+
+	double dist = geo.delta(&site->getGeo())
+				/ glass3::util::GlassMath::k_DegreesToRadians;
+
+	// log it
+	glass3::util::Logger::log(
+			"debug",
+			"CHypo::calculateResidual: Calculated residual: " + std::to_string(tRes)
+			+ "; Phase: " + phaseName + "; for Pick: " + pick->getID()
+			+ "; tObs: " + std::to_string(tObs) + "; dist: " + std::to_string(dist)
+			+ " and Hypo: " + getID());
+	*/
+
 	return (tRes);
 }
 
@@ -3077,7 +3100,7 @@ void CHypo::calculateStatistics() {
 	// pick distances and azimuths
 	std::vector<double> dis;
 	std::vector<double> azm;
-	m_iTeleseismicPhaseCount = 0;
+	int telePhaseCount = 0;
 	for (auto pick : m_vPickData) {
 		// get the site
 		std::shared_ptr<CSite> site = pick->getSite();
@@ -3088,7 +3111,7 @@ void CHypo::calculateStatistics() {
 
 		// count phases past teleseismic distance
 		if (delta >= CGlass::getTeleseismicDistanceLimit()) {
-			m_iTeleseismicPhaseCount++;
+			telePhaseCount++;
 		}
 
 		// add to distance vactor
@@ -3102,6 +3125,8 @@ void CHypo::calculateStatistics() {
 		azm.push_back(azimuth);
 	}
 
+	m_iTeleseismicPhaseCount = telePhaseCount;
+
 	// Calculate distance standard deviation. Note that the denominator is N
 	// and not N-1, since a mean of 0 is pre-ordained.
 	// The skewness is also 0, since the distribution
@@ -3146,6 +3171,17 @@ void CHypo::calculateStatistics() {
 
 	// compute gap
 	m_dGap = calculateGap(m_dLatitude, m_dLongitude, m_dDepth);
+
+	glass3::util::Logger::log("debug",
+		"CHypo::calculateStatistics: ID: " + getID()
+		+ "; m_dDistanceSD: " + std::to_string(m_dDistanceSD)
+		+ "; m_dMedianDistance: " + std::to_string(m_dMedianDistance)
+		+ "; m_dMinDistance: " + std::to_string(m_dMinDistance)
+		+ "; m_dAssociationDistanceCutoff: "
+		+ std::to_string(m_dAssociationDistanceCutoff)
+		+ "; m_iTeleseismicPhaseCount: "
+		+ std::to_string(m_iTeleseismicPhaseCount)
+		+ "; m_dGap: " + std::to_string(m_dGap));
 }
 
 // ---------------------------------------------------------trap

glasscore/glasslib/src/HypoList.cpp

@@ -24,7 +24,7 @@
 namespace glasscore {
 
 // constants
-const int CHypoList::k_nHypoSearchPastDurationForPick;
+constexpr double CHypoList::k_nHypoSearchPastDurationForPick;
 const int CHypoList::k_nMaxAllowableHypoCountDefault;
 const unsigned int CHypoList::k_nNumberOfMergeAnnealIterations;
 constexpr double CHypoList::k_dFinalMergeAnnealTimeStepSize;
@@ -122,7 +122,7 @@ bool CHypoList::associateData(std::shared_ptr<CPick> pk) {
 
 	// compute the list of hypos to associate with
 	// (a potential hypo must be before the pick we're associating)
-	// use the pick time minus 2400 seconds to compute the starting index
+	// use the pick time minus 3600 seconds to compute the starting index
 	// NOTE: Hard coded time delta
 	std::vector<std::weak_ptr<CHypo>> hypoList = getHypos(
 			pk->getTPick() - k_nHypoSearchPastDurationForPick, pk->getTPick());
@@ -473,7 +473,7 @@ bool CHypoList::processHypo(std::shared_ptr<CHypo> hyp) {
 			.count();
 
 	// Search for any associable picks that match hypo in the pick list
-	// NOTE: This uses the hard coded 2400 second scavenge duration default
+	// NOTE: This uses the hard coded 3600 second scavenge duration default
 	if (CGlass::getPickList()->scavenge(hyp)) {
 		// we should report this hypo since it has changed
 		breport = true;

glasscore/testdata/initialize.d

@@ -1 +1 @@
-{"Cmd": "Initialize","MaximumNumberOfPicks": 10000,"MaximumNumberOfPicksPerSite": 30,"MaximumNumberOfCorrelations": 1000,"MaximumNumberOfHypos": 250,"PickDuplicateWindow": 2.5,"NumberOfNucleationThreads": 5,"NumberOfHypoThreads": 3,"NumberOfWebThreads": 1,"SiteHoursWithoutPicking": 6,"SiteLookupInterval": 6,"SiteMaximumPicksPerHour": 200,"Params": {"NucleationStackThreshold": 0.5,"NucleationDataCountThreshold": 10,"AssociationStandardDeviationCutoff": 3.0,"PruningStandardDeviationCutoff": 3.0,"PickAffinityExponentialFactor": 2.5,"DistanceCutoffFactor": 5.0,"DistanceCutoffRatio": 0.8,"DistanceCutoffMinimum": 30.0,"HypoProcessCountLimit": 25,"CorrelationTimeWindow": 2.5,"CorrelationDistanceWindow": 0.5,"CorrelationCancelAge": 900,"BeamMatchingAzimuthWindow" : 22.5,"ReportingStackThreshold": 0.5,"ReportingDataThreshold":5,"EventFragmentDepthThreshold": 550.0,"EventFragmentAzimuthThreshold": 270.0},"PickClassification":{"NoiseClassificationThreshold": 0.75,"PhaseClassificationThreshold": 0.75,"DistanceClassificationThreshold": 0.75,"DistanceClassificationClasses": [0.0, 1.5, 10.0, 30.0],"DistanceClassificationClassesLowerBound": [0.0, 1.0, 5.0, 20.0],"DistanceClassificationClassesUpperBound": [4.0, 15.0, 35.0, 180.0],"AzimuthClassificationThreshold":0.75,"AzimuthClassificationUncertainty":90.0},"DefaultNucleationPhase": {"PhaseName": "P","TravFile": "./testdata/P.trv"},"AssociationPhases": [{"PhaseName": "P","Range": [ 0, 0, 120, 180 ],"TravFile": "./testdata/P.trv"},{"PhaseName": "S","Assoc": [ 10, 90 ],"TravFile": "./testdata/S.trv"}],"TestTravelTimes":false,"UseL1ResidualLocator":false,"PlottingInfo":{"graphicsOut":false,"graphicsStepKM":1.0,"graphicsSteps":100,"graphicsOutFolder":"./"}}
+{"Cmd": "Initialize","MaximumNumberOfPicks": 10000,"MaximumNumberOfPicksPerSite": 30,"MaximumNumberOfCorrelations": 1000,"MaximumNumberOfHypos": 250,"PickDuplicateWindow": 2.5,"NumberOfNucleationThreads": 5,"NumberOfHypoThreads": 3,"NumberOfWebThreads": 1,"SiteHoursWithoutPicking": 6,"SiteLookupInterval": 6,"SiteMaximumPicksPerHour": 200,"Params": {"NucleationStackThreshold": 0.5,"NucleationDataCountThreshold": 10,"AssociationStandardDeviationCutoff": 3.0,"PruningStandardDeviationCutoff": 3.0,"PickAffinityExponentialFactor": 2.5,"DistanceCutoffFactor": 5.0,"DistanceCutoffRatio": 0.8,"DistanceCutoffMinimum": 30.0,"HypoProcessCountLimit": 25,"CorrelationTimeWindow": 2.5,"CorrelationDistanceWindow": 0.5,"CorrelationCancelAge": 900,"BeamMatchingAzimuthWindow" : 22.5,"ReportingStackThreshold": 0.5,"ReportingDataThreshold":5,"EventFragmentDepthThreshold": 550.0,"EventFragmentAzimuthThreshold": 270.0},"PickClassification":{"NoiseClassificationThreshold": 0.75,"PhaseClassificationThreshold": 0.75,"DistanceClassificationThreshold": 0.75,"DistanceClassificationClasses": [0.0, 1.5, 10.0, 30.0],"DistanceClassificationClassesLowerBound": [0.0, 1.0, 5.0, 20.0],"DistanceClassificationClassesUpperBound": [4.0, 15.0, 35.0, 180.0],"AzimuthClassificationThreshold":0.75,"AzimuthClassificationUncertainty":90.0},"DefaultNucleationPhase": {"PhaseName": "P","TravFile": "./testdata/P.trv"},"AssociationPhases": [{"PhaseName": "P","Assoc": [ 0, 180 ],"TravFile": "./testdata/P.trv"},{"PhaseName": "S","Assoc": [ 10, 90 ],"TravFile": "./testdata/S.trv"}],"TestTravelTimes":false,"UseL1ResidualLocator":false,"PlottingInfo":{"graphicsOut":false,"graphicsStepKM":1.0,"graphicsSteps":100,"graphicsOutFolder":"./"}}

glasscore/tests/hypo_unittest.cpp

@@ -62,7 +62,7 @@
 #define ANNEAL_TIME 3648515732.3372693
 #define ANNEAL_BAYES 31.352108346757745
 
-#define LOCALIZE_LATITUDE 42.011810141430075
+#define LOCALIZE_LATITUDE 42.039847535832415
 #define LOCALIZE_LONGITUDE -119.36313949833216
 #define LOCALIZE_DEPTH 10.038057115005586
 #define LOCALIZE_TIME 3648515732.3254061

glasscore/traveltime/include/TTT.h

@@ -45,7 +45,7 @@ class CTTT {
 	 */
 	~CTTT();
 
-/**
+	/**
 	 * \brief Write out travel times to files
 	 * 
 	 * This function prints the available travel time points for the travel
@@ -175,7 +175,6 @@ class CTTT {
 	 */
 	// double testTravelTimes(std::string phase);
 
-
 	// constants
 	/**
 	 * \brief A value representing a travel time that is too large to be valid
@@ -205,7 +204,6 @@ class CTTT {
 	 */
 	bool m_bPublishable;
 
-
 	/**
 	 * \brief glass3::util::Geo object containing current
 	 * geographic location of source(as in source/receiver). Set by setOrigin()

glasscore/traveltime/src/TTT.cpp

@@ -231,11 +231,11 @@ double CTTT::T(double delta, std::string phase) {
 double CTTT::T(glass3::util::Geo *geo, double tObserved) {
 	// Find Phase with least residual, returns time
 
-	double bestTraveltime;
-	std::string bestPhase;
-	bool useForLocations;
-	bool publishable;
-	double minResidual = k_dTTTooLargeToBeValid;
+	double bestTraveltime = CTravelTime::k_dTravelTimeInvalid;
+	std::string bestPhase = "?";
+	bool useForLocations = true;
+	bool publishable = true;
+	double bestResidual = k_dTTTooLargeToBeValid;
 
 	// for each phase
 	for (int i = 0; i < m_iNumTravelTimes; i++) {
@@ -249,7 +249,7 @@ double CTTT::T(glass3::util::Geo *geo, double tObserved) {
 		double traveltime = aTrv->T(geo);
 
 		// check traveltime
-		if (traveltime < 0.0) {
+		if (traveltime <= CTravelTime::k_dTravelTimeInvalid) {
 			continue;
 		}
 
@@ -264,7 +264,7 @@ double CTTT::T(glass3::util::Geo *geo, double tObserved) {
 
 		// check to see if phase is associable
 		// based on maximum assoc distance, if present
-		if (m_adMaximumAssociationValues[i] > 0) {
+		if (m_adMaximumAssociationValues[i] >= 0) {
 			if (aTrv->m_dDelta > m_adMaximumAssociationValues[i]) {
 				// this phase is not associable  at this distance
 				continue;
@@ -275,10 +275,10 @@ double CTTT::T(glass3::util::Geo *geo, double tObserved) {
 		double residual = std::abs(tObserved - traveltime);
 
 		// check to see if this residual is better than the previous
-		//  best
-		if (residual < minResidual) {
+		// best
+		if (residual < bestResidual) {
 			// this is the new best travel time
-			minResidual = residual;
+			bestResidual = residual;
 			bestPhase = aTrv->m_sPhase;
 			bestTraveltime = traveltime;
 			useForLocations = aTrv->m_bUseForLocations;
@@ -287,7 +287,7 @@ double CTTT::T(glass3::util::Geo *geo, double tObserved) {
 	}
 
 	// check to see if minimum residual is valid
-	if (minResidual < k_dTTTooLargeToBeValid) {
+	if (bestResidual < k_dTTTooLargeToBeValid) {
 		m_sPhase = bestPhase;
 		m_bUseForLocations = useForLocations;
 		m_bPublishable = publishable;

neic-glass3.code-workspace

@@ -81,7 +81,12 @@
 			"complex": "cpp",
 			"*.tcc": "cpp",
 			"cinttypes": "cpp",
-			"hashtable": "cpp"
+			"hashtable": "cpp",
+			"__errc": "cpp",
+			"__node_handle": "cpp",
+			"bit": "cpp",
+			"condition_variable": "cpp",
+			"optional": "cpp"
 		}
 	}
 }