From 026772d971eaf5afc319c3495e88a1bde965c7b8 Mon Sep 17 00:00:00 2001 From: William L Yeck Date: Wed, 23 May 2018 11:47:28 -0600 Subject: [PATCH] Azimuthal Gap Downweighting option for local grids --- glasscore/glasslib/include/Hypo.h | 25 ++++++++- glasscore/glasslib/include/Web.h | 18 ++++++- glasscore/glasslib/src/Hypo.cpp | 79 ++++++++++++++++++++++++++--- glasscore/glasslib/src/HypoList.cpp | 2 +- glasscore/glasslib/src/Node.cpp | 9 ++-- glasscore/glasslib/src/Web.cpp | 43 +++++++++++++--- glasscore/tests/web_unittest.cpp | 3 +- 7 files changed, 156 insertions(+), 23 deletions(-) diff --git a/glasscore/glasslib/include/Hypo.h b/glasscore/glasslib/include/Hypo.h index d9aad97d..4297cd87 100644 --- a/glasscore/glasslib/include/Hypo.h +++ b/glasscore/glasslib/include/Hypo.h @@ -81,7 +81,8 @@ class CHypo { std::string web, double bayes, double thresh, int cut, std::shared_ptr firstTrav, std::shared_ptr secondTrav, - std::shared_ptr ttt, double resolution = 100); + std::shared_ptr ttt, double resolution = 100, + double aziTaper = 360.); /** * \brief CHypo alternate constructor @@ -148,7 +149,7 @@ class CHypo { std::shared_ptr firstTrav, std::shared_ptr secondTrav, std::shared_ptr ttt, double resolution = - 100); + 100, double aziTap = 360.); /** * \brief Add pick reference to this hypo @@ -478,6 +479,15 @@ class CHypo { double tStart, double tStop, int nucleate = 0); + /** + * Calculates azimuthal gap for a proposed location + * + * \param lat - latitude of test location + * \param lon - longitude of test location + * \param z - depth of test location + */ + double gap(double lat, double lon, double z); + /** * Gets the stack of associated arrivals at location * @@ -549,6 +559,12 @@ class CHypo { */ void trap(); + /** + * \brief aziTaper + * \return the aziTaper + */ + double getAziTaper() const; + /** * \brief Latitude getter * \return the latitude @@ -853,6 +869,11 @@ class CHypo { */ double searchVals[5]; + /** + * \brief where to taper bayesVal from azi gap + */ + double aziTaper; + /** * \brief An integer value containing this hypo's processing cycle count */ diff --git a/glasscore/glasslib/include/Web.h b/glasscore/glasslib/include/Web.h index 826b33c6..5fd3d9a7 100644 --- a/glasscore/glasslib/include/Web.h +++ b/glasscore/glasslib/include/Web.h @@ -92,7 +92,7 @@ class CWeb { std::shared_ptr firstTrav, std::shared_ptr secondTrav, int numThreads = 0, int sleepTime = 100, - int checkInterval = 60); + int checkInterval = 60, double aziTaper = 360.); /** * \brief CWeb destructor @@ -157,7 +157,8 @@ class CWeb { int numNucleate, int resolution, int numRows, int numCols, int numZ, bool update, std::shared_ptr firstTrav, - std::shared_ptr secondTrav); + std::shared_ptr secondTrav, + double aziTap); /** * \brief Generate a local detection grid @@ -338,6 +339,13 @@ class CWeb { */ void addJob(std::function newjob); + /** + * \brief aziTapre Getter + * \return double with the azi taper start + */ + double getAziTaper() const; + + /** * \brief CGlass getter * \return the CGlass pointer @@ -495,6 +503,12 @@ class CWeb { */ double dResolution; + /** + * \brief A double which describes where the locator should start + * down weighting for azimuthal gap + **/ + double aziTaper = 360.; + /** * \brief A boolean flag that stores whether to update this web when a * station has changed. diff --git a/glasscore/glasslib/src/Hypo.cpp b/glasscore/glasslib/src/Hypo.cpp index 35e77bf4..ca7d3d1e 100644 --- a/glasscore/glasslib/src/Hypo.cpp +++ b/glasscore/glasslib/src/Hypo.cpp @@ -52,13 +52,15 @@ CHypo::CHypo(double lat, double lon, double z, double time, std::string pid, std::string web, double bayes, double thresh, int cut, std::shared_ptr firstTrav, std::shared_ptr secondTrav, - std::shared_ptr ttt, double resolution) { + std::shared_ptr ttt, double resolution, + double aziTaper) { + // seed the random number generator std::random_device randomDevice; m_RandomGenerator.seed(randomDevice()); if (!initialize(lat, lon, z, time, pid, web, bayes, thresh, cut, firstTrav, - secondTrav, ttt, resolution)) { + secondTrav, ttt, resolution, aziTaper)) { clear(); } } @@ -93,7 +95,8 @@ CHypo::CHypo(std::shared_ptr trigger, trigger->getWeb()->getThresh(), trigger->getWeb()->getNucleate(), trigger->getWeb()->getTrv1(), trigger->getWeb()->getTrv2(), - ttt, trigger->getResolution())) { + ttt, trigger->getResolution(), + trigger->getWeb()->getAziTaper())) { clear(); } } @@ -412,11 +415,16 @@ void CHypo::annealingLocate(int nIter, double dStart, double dStop, return; } + // taper to lower val if large azimuthal gap + glassutil::CTaper taperGap; + taperGap = glassutil::CTaper(0., 0., aziTaper, 360.); + // these hold the values of the initial, current, and best stack location double valStart = 0; double valBest = 0; // calculate the value of the stack at the current location - valStart = getBayes(dLat, dLon, dZ, tOrg, nucleate); + valStart = getBayes(dLat, dLon, dZ, tOrg, nucleate) + * taperGap.Val(gap(dLat,dLon,dZ)); char sLog[1024]; @@ -483,7 +491,8 @@ void CHypo::annealingLocate(int nIter, double dStart, double dStop, double oT = tOrg + dt; // get the stack value for this hypocenter - double val = getBayes(xlat, xlon, xz, oT, nucleate); + double val = getBayes(xlat, xlon, xz, oT, nucleate) + * taperGap.Val(gap(xlat,xlon,xz)); // if testing locator print iteration if (pGlass->getTestLocator()) { @@ -497,6 +506,7 @@ void CHypo::annealingLocate(int nIter, double dStart, double dStop, // is this stacked bayesian value (val) better than the previous best // (valBest) + if (val > valBest || (val > dThresh && (valBest - val) @@ -1179,6 +1189,55 @@ std::shared_ptr CHypo::expire(bool send) { return (expire); } +// ---------------------------------------------------------Gap +double CHypo::gap(double lat, double lon, double z) { + // lock mutex for this scope + std::lock_guard guard(hypoMutex); + + // set up a geographic object for this hypo + glassutil::CGeo geo; + geo.setGeographic(lat, lon, EARTHRADIUSKM - z); + + // create and populate vectors containing the + // pick distances and azimuths + std::vector azm; + + for (auto pick : vPick) { + // get the site + std::shared_ptr site = pick->getSite(); + + // compute the azimuth + double azimuth = geo.azimuth(&site->getGeo()) / DEG2RAD; + + // add to azimuth vector + azm.push_back(azimuth); + } + + int nazm = azm.size(); + + if(nazm <= 1) + { + return 360.; + } + // sort the azimuths + sort(azm.begin(), azm.end()); + + // add the first (smallest) azimuth to the end by adding 360 + azm.push_back(azm.front() + 360.0); + + // compute gap + double tempGap = 0.0; + + for (int i = 0; i < nazm; i++) { + double gap = azm[i + 1] - azm[i]; + if (gap > tempGap) { + tempGap = gap; + } + } + + return tempGap; +} + // ---------------------------------------------------------Gaussian double CHypo::gauss(double avg, double std) { // generate Gaussian pseudo-random number using the @@ -1227,6 +1286,12 @@ double CHypo::getResidual(std::shared_ptr pick) { return tRes; } +double CHypo::getAziTaper() const { + std::lock_guard hypoGuard(hypoMutex); + return (aziTaper); +} + + double CHypo::getBayes() const { std::lock_guard hypoGuard(hypoMutex); return (dBayes); @@ -1903,7 +1968,7 @@ bool CHypo::initialize(double lat, double lon, double z, double time, std::shared_ptr firstTrav, std::shared_ptr secondTrav, std::shared_ptr ttt, - double resolution) { + double resolution, double aziTap) { // lock mutex for this scope std::lock_guard guard(hypoMutex); @@ -1917,7 +1982,7 @@ bool CHypo::initialize(double lat, double lon, double z, double time, sWebName = web; dBayes = bayes; dBayesInitial = bayes; - + aziTaper = aziTap; dThresh = thresh; nCut = cut; dRes = resolution; diff --git a/glasscore/glasslib/src/HypoList.cpp b/glasscore/glasslib/src/HypoList.cpp index 99a00c99..4f494b31 100644 --- a/glasscore/glasslib/src/HypoList.cpp +++ b/glasscore/glasslib/src/HypoList.cpp @@ -1133,7 +1133,7 @@ bool CHypoList::mergeCloseEvents(std::shared_ptr hypo) { glassutil::CPid::pid(), "Merged Hypo", 0.0, hypo->getThresh(), hypo->getCut(), hypo->getTrv1(), hypo->getTrv2(), pGlass->getTTT(), - distanceCut * DEG2KM); + distanceCut * DEG2KM, hypo->getAziTaper()); // set hypo glass pointer and such hypo3->setGlass(pGlass); diff --git a/glasscore/glasslib/src/Node.cpp b/glasscore/glasslib/src/Node.cpp index 466fdc6b..8937d8e6 100644 --- a/glasscore/glasslib/src/Node.cpp +++ b/glasscore/glasslib/src/Node.cpp @@ -509,11 +509,14 @@ std::string CNode::getSitesString() { for (const auto &link : vSite) { // get the site std::shared_ptr currentSite = std::get< LINK_PTR>(link); + double lat, lon, r; + + currentSite->getGeo().getGeographic(&lat, &lon, &r); siteString += sPid + "," + currentSite->getScnl() + "," - + std::to_string(currentSite->getGeo().dLat) + ";" - + std::to_string(currentSite->getGeo().dLon) + ";" - + std::to_string(currentSite->getGeo().dRad) + "\n"; + + std::to_string(lat) + ";" + + std::to_string(lon) + ";" + + std::to_string(r) + "\n"; } return (siteString); diff --git a/glasscore/glasslib/src/Web.cpp b/glasscore/glasslib/src/Web.cpp index 8cf6ac33..87180221 100644 --- a/glasscore/glasslib/src/Web.cpp +++ b/glasscore/glasslib/src/Web.cpp @@ -82,7 +82,7 @@ CWeb::CWeb(std::string name, double thresh, int numDetect, int numNucleate, int resolution, int numRows, int numCols, int numZ, bool update, std::shared_ptr firstTrav, std::shared_ptr secondTrav, int numThreads, - int sleepTime, int checkInterval) { + int sleepTime, int checkInterval, double aziTaper) { // setup threads if (numThreads > 0) { m_bRunProcessLoop = true; @@ -98,7 +98,7 @@ CWeb::CWeb(std::string name, double thresh, int numDetect, int numNucleate, clear(); initialize(name, thresh, numDetect, numNucleate, resolution, numRows, - numCols, numZ, update, firstTrav, secondTrav); + numCols, numZ, update, firstTrav, secondTrav, aziTaper); m_StatusMutex.lock(); m_ThreadStatusMap.clear(); @@ -197,7 +197,8 @@ bool CWeb::initialize(std::string name, double thresh, int numDetect, int numNucleate, int resolution, int numRows, int numCols, int numZ, bool update, std::shared_ptr firstTrav, - std::shared_ptr secondTrav) { + std::shared_ptr secondTrav, + double aziTap) { std::lock_guard webGuard(m_WebMutex); sName = name; @@ -211,7 +212,7 @@ bool CWeb::initialize(std::string name, double thresh, int numDetect, bUpdate = update; pTrv1 = firstTrav; pTrv2 = secondTrav; - + aziTaper = aziTap; // done return (true); } @@ -283,6 +284,7 @@ bool CWeb::global(std::shared_ptr com) { int zs = 0; bool saveGrid = false; bool update = false; + double aziTaper = 360.; // get grid configuration from json // name @@ -330,6 +332,13 @@ bool CWeb::global(std::shared_ptr com) { thresh = (*com)["Thresh"].ToDouble(); } + // sets the aziTaper value + if ((*com).HasKey("AzimuthGapTaper") + && ((*com)["AzimuthGapTaper"].GetType() + == json::ValueType::DoubleVal)) { + aziTaper = (*com)["AzimuthGapTaper"].ToDouble(); + } + // Node resolution for this Global grid if (((*com).HasKey("Resolution")) && ((*com)["Resolution"].GetType() == json::ValueType::DoubleVal)) { @@ -380,7 +389,7 @@ bool CWeb::global(std::shared_ptr com) { // init, note global doesn't use nRow or nCol initialize(name, thresh, detect, nucleate, resol, 0, 0, zs, update, pTrv1, - pTrv2); + pTrv2, aziTaper); // generate site and network filter lists genSiteFilters(com); @@ -528,6 +537,7 @@ bool CWeb::grid(std::shared_ptr com) { double lon = 0; int rows = 0; int cols = 0; + double aziTaper = 360.; std::vector zzz; int zs = 0; bool saveGrid = false; @@ -645,6 +655,13 @@ bool CWeb::grid(std::shared_ptr com) { return (false); } + // sets the aziTaper value + if ((*com).HasKey("AzimuthGapTaper") + && ((*com)["AzimuthGapTaper"].GetType() + == json::ValueType::DoubleVal)) { + aziTaper = (*com)["AzimuthGapTaper"].ToDouble(); + } + // whether to create a file detailing the node configuration for // this grid if ((*com).HasKey("SaveGrid")) { @@ -668,7 +685,7 @@ bool CWeb::grid(std::shared_ptr com) { // initialize initialize(name, thresh, detect, nucleate, resol, rows, cols, zs, update, - pTrv1, pTrv2); + pTrv1, pTrv2, aziTaper); // generate site and network filter lists genSiteFilters(com); @@ -830,6 +847,7 @@ bool CWeb::grid_explicit(std::shared_ptr com) { int nN = 0; bool saveGrid = false; bool update = false; + double aziTaper = 360.; std::vector> nodes; double resol = 0; @@ -880,6 +898,13 @@ bool CWeb::grid_explicit(std::shared_ptr com) { thresh = (*com)["Thresh"].ToDouble(); } + // sets the aziTaper value + if ((*com).HasKey("AzimuthGapTaper") + && ((*com)["AzimuthGapTaper"].GetType() + == json::ValueType::DoubleVal)) { + aziTaper = (*com)["AzimuthGapTaper"].ToDouble(); + } + // Node resolution for this grid if (((*com).HasKey("Resolution")) && ((*com)["Resolution"].GetType() == json::ValueType::DoubleVal)) { @@ -958,7 +983,7 @@ bool CWeb::grid_explicit(std::shared_ptr com) { // initialize initialize(name, thresh, detect, nucleate, resol, 0., 0., 0., update, pTrv1, - pTrv2); + pTrv2, aziTaper); // generate site and network filter lists genSiteFilters(com); @@ -1979,6 +2004,10 @@ bool CWeb::hasSite(std::shared_ptr site) { return (false); } +double CWeb::getAziTaper() const { + return (aziTaper); +} + const CSiteList* CWeb::getSiteList() const { std::lock_guard webGuard(m_WebMutex); return (pSiteList); diff --git a/glasscore/tests/web_unittest.cpp b/glasscore/tests/web_unittest.cpp index acfda969..18496eec 100644 --- a/glasscore/tests/web_unittest.cpp +++ b/glasscore/tests/web_unittest.cpp @@ -34,6 +34,7 @@ #define NUMDETECT 5 #define NUMNUCLEATE 4 #define RESOLUTION 100.0 +#define AZIGAP 360. #define NUMROWS 3 #define NUMCOLS 4 #define NUMZ 1 @@ -205,7 +206,7 @@ TEST(WebTest, Initialize) { NUMCOLS, NUMZ, UPDATE, - nullTrav, nullTrav); + nullTrav, nullTrav, AZIGAP); // name ASSERT_STREQ(std::string(NAME).c_str(), testWeb.getName().c_str())<<