diff --git a/src/base/Makefile b/src/base/Makefile
index a4ca203..891f1d0 100644
--- a/src/base/Makefile
+++ b/src/base/Makefile
@@ -33,7 +33,8 @@ FILES= Announce.cpp \
PolynomialInterp.cpp \
LegendrePolynomial.cpp \
FunctionTimer.cpp \
- ShpFile.cpp
+ ShpFile.cpp \
+ ThresholdOp.cpp
LIB_TARGET= libextremesbase.a
diff --git a/src/base/ThresholdOp.cpp b/src/base/ThresholdOp.cpp
new file mode 100644
index 0000000..49674dc
--- /dev/null
+++ b/src/base/ThresholdOp.cpp
@@ -0,0 +1,612 @@
+///////////////////////////////////////////////////////////////////////////////
+///
+/// \file ThresholdOp.cpp
+/// \author Paul Ullrich
+/// \version February 9, 2024
+///
+///
+/// Copyright 2024 Paul Ullrich
+///
+/// This file is distributed as part of the Tempest source code package.
+/// Permission is granted to use, copy, modify and distribute this
+/// source code and its documentation under the terms of the GNU General
+/// Public License. This software is provided "as is" without express
+/// or implied warranty.
+///
+
+#include "ThresholdOp.h"
+#include "SimpleGrid.h"
+#include "DataArray1D.h"
+
+#include
+#include
+#include
+
+///////////////////////////////////////////////////////////////////////////////
+// ThresholdOp
+///////////////////////////////////////////////////////////////////////////////
+
+void ThresholdOp::Parse(
+ VariableRegistry & varreg,
+ const std::string & strOp
+) {
+ // Read mode
+ enum {
+ ReadMode_Op,
+ ReadMode_Value,
+ ReadMode_Distance,
+ ReadMode_Invalid
+ } eReadMode = ReadMode_Op;
+
+ // Parse variable
+ int iLast = varreg.FindOrRegisterSubStr(strOp, &m_varix) + 1;
+
+ // Loop through string
+ for (int i = iLast; i <= strOp.length(); i++) {
+
+ // Comma-delineated
+ if ((i == strOp.length()) || (strOp[i] == ',')) {
+
+ std::string strSubStr =
+ strOp.substr(iLast, i - iLast);
+
+ // Read in operation
+ if (eReadMode == ReadMode_Op) {
+ if (strSubStr == ">") {
+ m_eOp = GreaterThan;
+ } else if (strSubStr == "<") {
+ m_eOp = LessThan;
+ } else if (strSubStr == ">=") {
+ m_eOp = GreaterThanEqualTo;
+ } else if (strSubStr == "<=") {
+ m_eOp = LessThanEqualTo;
+ } else if (strSubStr == "=") {
+ m_eOp = EqualTo;
+ } else if (strSubStr == "!=") {
+ m_eOp = NotEqualTo;
+ } else {
+ _EXCEPTION1("Threshold invalid operation \"%s\"",
+ strSubStr.c_str());
+ }
+
+ iLast = i + 1;
+ eReadMode = ReadMode_Value;
+
+ // Read in value
+ } else if (eReadMode == ReadMode_Value) {
+ m_dThresholdValue = atof(strSubStr.c_str());
+
+ iLast = i + 1;
+ eReadMode = ReadMode_Distance;
+
+ // Read in distance to point that satisfies threshold
+ } else if (eReadMode == ReadMode_Distance) {
+ m_dDistanceDeg = atof(strSubStr.c_str());
+
+ iLast = i + 1;
+ eReadMode = ReadMode_Invalid;
+
+ // Invalid
+ } else if (eReadMode == ReadMode_Invalid) {
+ _EXCEPTION1("\nInsufficient entries in threshold op \"%s\""
+ "\nRequired: \","
+ ",,\"",
+ strOp.c_str());
+ }
+ }
+ }
+
+ if (eReadMode != ReadMode_Invalid) {
+ _EXCEPTION1("\nInsufficient entries in threshold op \"%s\""
+ "\nRequired: \",,,\"",
+ strOp.c_str());
+ }
+
+ if (m_dDistanceDeg < 0.0) {
+ _EXCEPTIONT("For threshold op, distance must be nonnegative");
+ }
+ if (m_dDistanceDeg > 180.0) {
+ _EXCEPTIONT("For threshold op, distance must be less than 180 degrees");
+ }
+
+ // Output announcement
+ std::string strDescription = varreg.GetVariableString(m_varix);
+ if (m_eOp == GreaterThan) {
+ strDescription += " is greater than ";
+ } else if (m_eOp == LessThan) {
+ strDescription += " is less than ";
+ } else if (m_eOp == GreaterThanEqualTo) {
+ strDescription += " is greater than or equal to ";
+ } else if (m_eOp == LessThanEqualTo) {
+ strDescription += " is less than or equal to ";
+ } else if (m_eOp == EqualTo) {
+ strDescription += " is equal to ";
+ } else if (m_eOp == NotEqualTo) {
+ strDescription += " is not equal to ";
+ }
+
+ char szBuffer[128];
+ snprintf(szBuffer, 128, "%f within %f degrees",
+ m_dThresholdValue, m_dDistanceDeg);
+ strDescription += szBuffer;
+
+ Announce("%s", strDescription.c_str());
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+template
+bool ThresholdOp::IsSatisfied(
+ const SimpleGrid & grid,
+ const DataArray1D & dataState,
+ const int ix0
+) const {
+ // Special case if m_dDistanceDeg is zero
+ if (m_dDistanceDeg < 1.0e-12) {
+ double dValue = dataState[ix0];
+
+ if (m_eOp == ThresholdOp::GreaterThan) {
+ if (dValue > m_dThresholdValue) {
+ return true;
+ }
+
+ } else if (m_eOp == ThresholdOp::LessThan) {
+ if (dValue < m_dThresholdValue) {
+ return true;
+ }
+
+ } else if (m_eOp == ThresholdOp::GreaterThanEqualTo) {
+ if (dValue >= m_dThresholdValue) {
+ return true;
+ }
+
+ } else if (m_eOp == ThresholdOp::LessThanEqualTo) {
+ if (dValue <= m_dThresholdValue) {
+ return true;
+ }
+
+ } else if (m_eOp == ThresholdOp::EqualTo) {
+ if (dValue == m_dThresholdValue) {
+ return true;
+ }
+
+ } else if (m_eOp == ThresholdOp::NotEqualTo) {
+ if (dValue != m_dThresholdValue) {
+ return true;
+ }
+
+ } else {
+ _EXCEPTIONT("Invalid operation");
+ }
+
+ return false;
+ }
+
+ // Queue of nodes that remain to be visited
+ std::queue queueNodes;
+ queueNodes.push(ix0);
+
+ // Set of nodes that have already been visited
+ std::set setNodesVisited;
+
+ // Latitude and longitude at the origin
+ double dLat0 = grid.m_dLat[ix0];
+ double dLon0 = grid.m_dLon[ix0];
+
+ // Loop through all latlon elements
+ while (queueNodes.size() != 0) {
+ int ix = queueNodes.front();
+ queueNodes.pop();
+
+ if (setNodesVisited.find(ix) != setNodesVisited.end()) {
+ continue;
+ }
+
+ setNodesVisited.insert(ix);
+
+ // Great circle distance to this element
+ double dLatThis = grid.m_dLat[ix];
+ double dLonThis = grid.m_dLon[ix];
+
+ double dR =
+ sin(dLat0) * sin(dLatThis)
+ + cos(dLat0) * cos(dLatThis) * cos(dLonThis - dLon0);
+
+ if (dR >= 1.0) {
+ dR = 0.0;
+ } else if (dR <= -1.0) {
+ dR = 180.0;
+ } else {
+ dR = 180.0 / M_PI * acos(dR);
+ }
+ if (dR != dR) {
+ _EXCEPTIONT("NaN value detected");
+ }
+
+ if ((ix != ix0) && (dR > m_dDistanceDeg)) {
+ continue;
+ }
+
+ // Value at this location
+ double dValue = dataState[ix];
+
+ // Apply operator
+ if (m_eOp == ThresholdOp::GreaterThan) {
+ if (dValue > m_dThresholdValue) {
+ return true;
+ }
+
+ } else if (m_eOp == ThresholdOp::LessThan) {
+ if (dValue < m_dThresholdValue) {
+ return true;
+ }
+
+ } else if (m_eOp == ThresholdOp::GreaterThanEqualTo) {
+ if (dValue >= m_dThresholdValue) {
+ return true;
+ }
+
+ } else if (m_eOp == ThresholdOp::LessThanEqualTo) {
+ if (dValue <= m_dThresholdValue) {
+ return true;
+ }
+
+ } else if (m_eOp == ThresholdOp::EqualTo) {
+ if (dValue == m_dThresholdValue) {
+ return true;
+ }
+
+ } else if (m_eOp == ThresholdOp::NotEqualTo) {
+ if (dValue != m_dThresholdValue) {
+ return true;
+ }
+
+ } else {
+ _EXCEPTIONT("Invalid operation");
+ }
+
+ // Add all neighbors of this point
+ for (int n = 0; n < grid.m_vecConnectivity[ix].size(); n++) {
+ queueNodes.push(grid.m_vecConnectivity[ix][n]);
+ }
+ }
+
+ return false;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// ThresholdOpTreeNode
+///////////////////////////////////////////////////////////////////////////////
+
+void ThresholdOpTreeNode::Parse(
+ VariableRegistry & varreg,
+ std::string strExpression,
+ bool fAllowZeroLengthExpression
+) {
+ if (m_vecChildren.size() != 0) {
+ _EXCEPTIONT("ThresholdOpTreeNode already initialized");
+ }
+
+ // Set to a leaf node by default
+ m_eType = AlwaysTrueNode;
+
+ // Pre-process expression
+ for (;;) {
+ // Remove whitespace on either side of expression
+ STLStringHelper::RemoveWhitespaceInPlace(strExpression);
+
+ // Zero length expressions are always true
+ if (strExpression.length() == 0) {
+ if (fAllowZeroLengthExpression) {
+ return;
+ } else {
+ _EXCEPTIONT("Zero-length sub-expression in threshold operator");
+ }
+ }
+
+ // Remove parentheses if they wrap whole expression, then re-process
+ if ((strExpression[0] == '(') && (strExpression[strExpression.length()-1] == ')')) {
+ int nParentheses = 0;
+ for (int i = 0; i < strExpression.length(); i++) {
+ if (strExpression[i] == '(') {
+ nParentheses++;
+ }
+ if (strExpression[i] == ')') {
+ if (nParentheses == 0) {
+ _EXCEPTION1("Malformed sub-expression in threshold operator \"%s\"", strExpression.c_str());
+ }
+ nParentheses--;
+ if ((nParentheses == 0) && (i != strExpression.length()-1)) {
+ nParentheses = -1;
+ break;
+ }
+ }
+ }
+ if (nParentheses == 0) {
+ strExpression = strExpression.substr(1,strExpression.length()-2);
+ continue;
+ }
+ }
+
+ break;
+ }
+
+ // Check for operators and balanced parentheses
+ // If there are any "or" operators in strExpression then vecOpPositionIx only contains "or" positions
+ // This is because "and" operations have precedence and so are evaluated in the children
+ std::vector vecOpPositionIx;
+ int nParentheses = 0;
+ for (int i = 0; i < strExpression.length(); i++) {
+ if (strExpression[i] == '(') {
+ nParentheses++;
+ }
+ if (strExpression[i] == ')') {
+ if (nParentheses == 0) {
+ _EXCEPTION1("Malformed sub-expression in threshold operator \"%s\"", strExpression.c_str());
+ }
+ nParentheses--;
+ }
+ if (nParentheses == 0) {
+ if (strExpression[i] == '|') {
+ if (m_eType != OrNode) {
+ m_eType = OrNode;
+ vecOpPositionIx.clear();
+ }
+ vecOpPositionIx.push_back(i);
+ }
+ if ((strExpression[i] == '&') || (strExpression[i] == ';')) {
+ if (m_eType != OrNode) {
+ m_eType = AndNode;
+ vecOpPositionIx.push_back(i);
+ }
+ }
+ }
+ }
+ if (nParentheses != 0) {
+ _EXCEPTION1("Malformed sub-expression in threshold operator \"%s\"", strExpression.c_str());
+ }
+
+ // Operator node; split operations
+ if (vecOpPositionIx.size() != 0) {
+ int iLast = 0;
+ for (int i = 0; i < vecOpPositionIx.size(); i++) {
+ _ASSERT(vecOpPositionIx[i] - iLast - 1 >= 0);
+
+ ThresholdOpTreeNode * pChild = new ThresholdOpTreeNode();
+ pChild->Parse(
+ varreg,
+ strExpression.substr(
+ iLast, vecOpPositionIx[i] - iLast),
+ false);
+ m_vecChildren.push_back(pChild);
+
+ iLast = vecOpPositionIx[i]+1;
+ }
+
+ ThresholdOpTreeNode * pChild = new ThresholdOpTreeNode();
+ pChild->Parse(
+ varreg,
+ strExpression.substr(iLast),
+ false);
+ m_vecChildren.push_back(pChild);
+
+ // Leaf node; evaluate as threshold operator
+ } else {
+ m_eType = ThresholdNode;
+ m_threshop.Parse(varreg, strExpression);
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+template
+bool ThresholdOpTreeNode::IsSatisfied(
+ VariableRegistry & varreg,
+ NcFileVector & vecFiles,
+ const SimpleGrid & grid,
+ const int ix0,
+ std::vector * pvecRejectedCount
+) const {
+ if (pvecRejectedCount != NULL) {
+ _ASSERT(pvecRejectedCount->size() == size());
+ }
+ if (m_eType == AlwaysTrueNode) {
+ return true;
+ }
+ if (m_eType == ThresholdNode) {
+
+ Variable & var = varreg.Get(m_threshop.m_varix);
+ var.LoadGridData(varreg, vecFiles, grid);
+ const DataArray1D & dataState = var.GetData();
+
+ bool fCurrentSatisfied = m_threshop.IsSatisfied(grid, dataState, ix0);
+ if (pvecRejectedCount != NULL) {
+ (*pvecRejectedCount)[0] += (fCurrentSatisfied)?(0):(1);
+ }
+ return fCurrentSatisfied;
+ }
+ if (m_eType == OrNode) {
+ _ASSERT(m_vecChildren.size() != 0);
+ bool fAnySatisfied = false;
+ for (int i = 0; i < m_vecChildren.size(); i++) {
+ bool fCurrentSatisfied = m_vecChildren[i]->IsSatisfied(varreg, vecFiles, grid, ix0);
+ if (pvecRejectedCount != NULL) {
+ (*pvecRejectedCount)[i] += (fCurrentSatisfied)?(0):(1);
+ }
+ fAnySatisfied = fAnySatisfied || fCurrentSatisfied;
+ }
+ return fAnySatisfied;
+ }
+ if (m_eType == AndNode) {
+ _ASSERT(m_vecChildren.size() != 0);
+ bool fAllSatisfied = true;
+ for (int i = 0; i < m_vecChildren.size(); i++) {
+ bool fCurrentSatisfied = m_vecChildren[i]->IsSatisfied(varreg, vecFiles, grid, ix0);
+ if (pvecRejectedCount != NULL) {
+ (*pvecRejectedCount)[i] += (fCurrentSatisfied)?(0):(1);
+ }
+ fAllSatisfied = fAllSatisfied && fCurrentSatisfied;
+ }
+ return fAllSatisfied;
+ }
+ _EXCEPTION();
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+template
+void ThresholdOpTreeNode::IsSatisfied(
+ VariableRegistry & varreg,
+ NcFileVector & vecFiles,
+ const SimpleGrid & grid,
+ DataArray1D & bTag
+) const {
+ _ASSERT(bTag.GetRows() == grid.GetSize());
+
+ if (m_eType == AlwaysTrueNode) {
+ return;
+ }
+ if (m_eType == ThresholdNode) {
+
+ Variable & var = varreg.Get(m_threshop.m_varix);
+ var.LoadGridData(varreg, vecFiles, grid);
+ const DataArray1D & dataState = var.GetData();
+
+ for (int ix = 0; ix < bTag.GetRows(); ix++) {
+ if (bTag[ix] == 0) {
+ continue;
+ }
+ if (!m_threshop.IsSatisfied(grid, dataState, ix)) {
+ bTag[ix] = 0;
+ }
+ }
+ return;
+ }
+ if (m_eType == OrNode) {
+ _ASSERT(m_vecChildren.size() != 0);
+ DataArray1D bTagCopy(bTag.GetRows());
+ for (int i = 0; i < m_vecChildren.size(); i++) {
+ std::memcpy(&(bTagCopy[0]), &(bTag[0]), bTag.GetRows() * sizeof(int));
+ for (int j = 0; j < bTagCopy.GetRows(); j++) {
+ if (bTagCopy[j] == (-1)) {
+ bTagCopy[j] = 1;
+ }
+ }
+ m_vecChildren[i]->IsSatisfied(varreg, vecFiles, grid, bTagCopy);
+ for (int j = 0; j < bTag.GetRows(); j++) {
+ if (bTagCopy[j] != 0) {
+ bTag[j] = (-1);
+ }
+ }
+ }
+ for (int j = 0; j < bTag.GetRows(); j++) {
+ if (bTag[j] == 1) {
+ bTag[j] = 0;
+ }
+ if (bTag[j] == (-1)) {
+ bTag[j] = 1;
+ }
+ }
+ return;
+ }
+ if (m_eType == AndNode) {
+ _ASSERT(m_vecChildren.size() != 0);
+ for (int i = 0; i < m_vecChildren.size(); i++) {
+ m_vecChildren[i]->IsSatisfied(varreg, vecFiles, grid, bTag);
+ }
+ return;
+ }
+ _EXCEPTION();
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+template
+void ThresholdOpTreeNode::IsSatisfied(
+ VariableRegistry & varreg,
+ NcFileVector & vecFiles,
+ const SimpleGrid & grid,
+ const std::set & setCandidates,
+ std::set & setNewCandidates
+) const {
+ if (m_eType == AlwaysTrueNode) {
+ setNewCandidates = setCandidates;
+ return;
+ }
+ if (m_eType == ThresholdNode) {
+ Variable & var = varreg.Get(m_threshop.m_varix);
+ var.LoadGridData(varreg, vecFiles, grid);
+ const DataArray1D & dataState = var.GetData();
+
+ for (auto it = setCandidates.begin(); it != setCandidates.end(); it++) {
+ if (m_threshop.IsSatisfied(grid, dataState, *it)) {
+ setNewCandidates.insert(*it);
+ }
+ }
+ _ASSERT(setNewCandidates.size() <= setCandidates.size());
+ return;
+ }
+ if (m_eType == OrNode) {
+ _ASSERT(m_vecChildren.size() != 0);
+ for (int i = 0; i < m_vecChildren.size(); i++) {
+ std::set setNewCandidatesChild;
+
+ m_vecChildren[i]->IsSatisfied(varreg, vecFiles, grid, setCandidates, setNewCandidatesChild);
+
+ for (auto it = setNewCandidatesChild.begin(); it != setNewCandidatesChild.end(); it++) {
+ setNewCandidates.insert(*it);
+ }
+ }
+ _ASSERT(setNewCandidates.size() <= setCandidates.size());
+ return;
+ }
+ if (m_eType == AndNode) {
+ _ASSERT(m_vecChildren.size() != 0);
+ m_vecChildren[0]->IsSatisfied(varreg, vecFiles, grid, setCandidates, setNewCandidates);
+ for (int i = 1; i < m_vecChildren.size(); i++) {
+ std::set setNewCandidatesChild;
+ for (auto it = setCandidates.begin(); it != setCandidates.end(); it++) {
+ m_vecChildren[i]->IsSatisfied(varreg, vecFiles, grid, setNewCandidates, setNewCandidatesChild);
+ }
+ setNewCandidates = setNewCandidatesChild;
+ }
+ _ASSERT(setNewCandidates.size() <= setCandidates.size());
+ return;
+ }
+ _EXCEPTION();
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Explicit template instantiation
+
+template bool ThresholdOp::IsSatisfied(
+ const SimpleGrid & grid,
+ const DataArray1D & dataState,
+ const int ix0
+) const;
+
+template bool ThresholdOpTreeNode::IsSatisfied(
+ VariableRegistry & varreg,
+ NcFileVector & vecFiles,
+ const SimpleGrid & grid,
+ const int ix0,
+ std::vector * pvecRejectedCount
+) const;
+
+template void ThresholdOpTreeNode::IsSatisfied(
+ VariableRegistry & varreg,
+ NcFileVector & vecFiles,
+ const SimpleGrid & grid,
+ DataArray1D & bTag
+) const;
+
+template void ThresholdOpTreeNode::IsSatisfied(
+ VariableRegistry & varreg,
+ NcFileVector & vecFiles,
+ const SimpleGrid & grid,
+ const std::set & setCandidates,
+ std::set & setNewCandidates
+) const;
+
+///////////////////////////////////////////////////////////////////////////////
+
diff --git a/src/base/ThresholdOp.h b/src/base/ThresholdOp.h
index b546d7f..b4f0212 100644
--- a/src/base/ThresholdOp.h
+++ b/src/base/ThresholdOp.h
@@ -19,8 +19,12 @@
#include "Announce.h"
#include "Variable.h"
+#include "NcFileVector.h"
+#include "SimpleGrid.h"
+#include "DataArray1D.h"
#include
+#include
///////////////////////////////////////////////////////////////////////////////
@@ -50,8 +54,8 @@ class ThresholdOp {
ThresholdOp() :
m_varix(InvalidVariableIndex),
m_eOp(GreaterThan),
- m_dValue(0.0),
- m_dDistance(0.0)
+ m_dThresholdValue(0.0),
+ m_dDistanceDeg(0.0)
{ }
public:
@@ -61,106 +65,18 @@ class ThresholdOp {
void Parse(
VariableRegistry & varreg,
const std::string & strOp
- ) {
- // Read mode
- enum {
- ReadMode_Op,
- ReadMode_Value,
- ReadMode_Distance,
- ReadMode_Invalid
- } eReadMode = ReadMode_Op;
-
- // Parse variable
- int iLast = varreg.FindOrRegisterSubStr(strOp, &m_varix) + 1;
-
- // Loop through string
- for (int i = iLast; i <= strOp.length(); i++) {
-
- // Comma-delineated
- if ((i == strOp.length()) || (strOp[i] == ',')) {
-
- std::string strSubStr =
- strOp.substr(iLast, i - iLast);
-
- // Read in operation
- if (eReadMode == ReadMode_Op) {
- if (strSubStr == ">") {
- m_eOp = GreaterThan;
- } else if (strSubStr == "<") {
- m_eOp = LessThan;
- } else if (strSubStr == ">=") {
- m_eOp = GreaterThanEqualTo;
- } else if (strSubStr == "<=") {
- m_eOp = LessThanEqualTo;
- } else if (strSubStr == "=") {
- m_eOp = EqualTo;
- } else if (strSubStr == "!=") {
- m_eOp = NotEqualTo;
- } else {
- _EXCEPTION1("Threshold invalid operation \"%s\"",
- strSubStr.c_str());
- }
-
- iLast = i + 1;
- eReadMode = ReadMode_Value;
-
- // Read in value
- } else if (eReadMode == ReadMode_Value) {
- m_dValue = atof(strSubStr.c_str());
-
- iLast = i + 1;
- eReadMode = ReadMode_Distance;
-
- // Read in distance to point that satisfies threshold
- } else if (eReadMode == ReadMode_Distance) {
- m_dDistance = atof(strSubStr.c_str());
-
- iLast = i + 1;
- eReadMode = ReadMode_Invalid;
-
- // Invalid
- } else if (eReadMode == ReadMode_Invalid) {
- _EXCEPTION1("\nInsufficient entries in threshold op \"%s\""
- "\nRequired: \","
- ",,\"",
- strOp.c_str());
- }
- }
- }
-
- if (eReadMode != ReadMode_Invalid) {
- _EXCEPTION1("\nInsufficient entries in threshold op \"%s\""
- "\nRequired: \",,,\"",
- strOp.c_str());
- }
-
- if (m_dDistance < 0.0) {
- _EXCEPTIONT("For threshold op, distance must be nonnegative");
- }
+ );
- // Output announcement
- std::string strDescription = varreg.GetVariableString(m_varix);
- if (m_eOp == GreaterThan) {
- strDescription += " is greater than ";
- } else if (m_eOp == LessThan) {
- strDescription += " is less than ";
- } else if (m_eOp == GreaterThanEqualTo) {
- strDescription += " is greater than or equal to ";
- } else if (m_eOp == LessThanEqualTo) {
- strDescription += " is less than or equal to ";
- } else if (m_eOp == EqualTo) {
- strDescription += " is equal to ";
- } else if (m_eOp == NotEqualTo) {
- strDescription += " is not equal to ";
- }
-
- char szBuffer[128];
- snprintf(szBuffer, 128, "%f within %f degrees",
- m_dValue, m_dDistance);
- strDescription += szBuffer;
-
- Announce("%s", strDescription.c_str());
- }
+public:
+ ///
+ /// Returns true if this threshold is satisfied, or false otherwise.
+ ///
+ template
+ bool IsSatisfied(
+ const SimpleGrid & grid,
+ const DataArray1D & dataState,
+ const int ix0
+ ) const;
public:
///
@@ -176,12 +92,146 @@ class ThresholdOp {
///
/// Threshold value.
///
- double m_dValue;
+ double m_dThresholdValue;
+
+ ///
+ /// Distance to search for threshold value (in degrees).
+ ///
+ double m_dDistanceDeg;
+};
+
+///////////////////////////////////////////////////////////////////////////////
+
+///
+/// A class for holding a tree node from a threshold operator evaluation tree.
+///
+class ThresholdOpTreeNode {
+
+public:
+ ///
+ /// Types of threshold operators.
+ ///
+ enum NodeType {
+ AlwaysTrueNode,
+ ThresholdNode,
+ AndNode,
+ OrNode
+ };
+
+public:
+ ///
+ /// Constructor.
+ ///
+ ThresholdOpTreeNode(
+ NodeType eType = AlwaysTrueNode
+ ) :
+ m_eType(eType)
+ { }
+
+ ///
+ /// Destructor.
+ ///
+ ~ThresholdOpTreeNode() {
+ for (auto pchild : m_vecChildren) {
+ delete pchild;
+ }
+ }
+
+ ///
+ /// Parse an expression and populate this node.
+ ///
+ void Parse(
+ VariableRegistry & varreg,
+ std::string strExpression,
+ bool fAllowZeroLengthExpression = true
+ );
+
+ ///
+ /// Returns true if this threshold is satisfied.
+ ///
+ template
+ bool IsSatisfied(
+ VariableRegistry & varreg,
+ NcFileVector & vecFiles,
+ const SimpleGrid & grid,
+ const int ix0,
+ std::vector * pvecRejectedCount = NULL
+ ) const;
+
+ ///
+ /// Performs IsSatisfied in place over all values in a binary mask.
+ ///
+ template
+ void IsSatisfied(
+ VariableRegistry & varreg,
+ NcFileVector & vecFiles,
+ const SimpleGrid & grid,
+ DataArray1D & bTag
+ ) const;
+
+ ///
+ /// Performs IsSatisfied over a set of candidate nodes, returning
+ /// the set of candidates that satisfy the criteria.
+ ///
+ template
+ void IsSatisfied(
+ VariableRegistry & varreg,
+ NcFileVector & vecFiles,
+ const SimpleGrid & grid,
+ const std::set & setCandidates,
+ std::set & setNewCandidates
+ ) const;
+
+ ///
+ /// If AlwaysTrueNode return 0.
+ /// If ThresholdNode return 1.
+ /// If AndNode or OrNode return the number of children.
+ ///
+ size_t size() const {
+ if (m_eType == AlwaysTrueNode) {
+ return 0;
+ }
+ if (m_eType == ThresholdNode) {
+ return 1;
+ }
+ return m_vecChildren.size();
+ }
+
+ ///
+ /// Get the name of this ThresholdOpTreeNode.
+ ///
+ std::string ToString(
+ VariableRegistry & varreg,
+ int ix
+ ) const {
+ _ASSERT((ix >= 0) && (ix < m_vecChildren.size()));
+ if (m_eType == AlwaysTrueNode) {
+ return 0;
+ }
+ if (m_eType == ThresholdNode) {
+ return varreg.GetVariableString(m_threshop.m_varix);
+ }
+ if (m_vecChildren[ix]->m_eType == ThresholdNode) {
+ return m_vecChildren[ix]->ToString(varreg, 0);
+ }
+ return std::string("Compound Expr. ") + std::to_string(ix);
+ }
+
+public:
+ ///
+ /// Node type.
+ ///
+ NodeType m_eType;
+
+ ///
+ /// Threshold operator at this node.
+ ///
+ ThresholdOp m_threshop;
///
- /// Distance to search for threshold value
+ /// Children nodes.
///
- double m_dDistance;
+ std::vector m_vecChildren;
};
///////////////////////////////////////////////////////////////////////////////
diff --git a/src/blobs/DetectBlobs.cpp b/src/blobs/DetectBlobs.cpp
index c7b7c73..b2dde41 100644
--- a/src/blobs/DetectBlobs.cpp
+++ b/src/blobs/DetectBlobs.cpp
@@ -45,162 +45,6 @@ typedef std::pair Point;
///////////////////////////////////////////////////////////////////////////////
-///
-/// Determine if the given field satisfies the threshold.
-///
-template
-bool SatisfiesThreshold(
- const SimpleGrid & grid,
- const DataArray1D & dataState,
- const int ix0,
- const ThresholdOp::Operation op,
- const double dTargetValue,
- const double dMaxDist
-) {
- // Special case if dMaxDist is zero
- if (dMaxDist < 1.0e-12) {
- double dValue = dataState[ix0];
-
- if (op == ThresholdOp::GreaterThan) {
- if (dValue > dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::LessThan) {
- if (dValue < dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::GreaterThanEqualTo) {
- if (dValue >= dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::LessThanEqualTo) {
- if (dValue <= dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::EqualTo) {
- if (dValue == dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::NotEqualTo) {
- if (dValue != dTargetValue) {
- return true;
- }
-
- } else {
- _EXCEPTIONT("Invalid operation");
- }
-
- }
-
- // Verify that dMaxDist is less than 180.0
- if (dMaxDist > 180.0) {
- _EXCEPTIONT("MaxDist must be less than 180.0");
- }
-
- // Queue of nodes that remain to be visited
- std::queue queueNodes;
- queueNodes.push(ix0);
-
- // Set of nodes that have already been visited
- std::set setNodesVisited;
-
- // Latitude and longitude at the origin
- double dLat0 = grid.m_dLat[ix0];
- double dLon0 = grid.m_dLon[ix0];
-
- // Loop through all latlon elements
- while (queueNodes.size() != 0) {
- int ix = queueNodes.front();
- queueNodes.pop();
-
- if (setNodesVisited.find(ix) != setNodesVisited.end()) {
- continue;
- }
-
- setNodesVisited.insert(ix);
-
- // Great circle distance to this element
- double dLatThis = grid.m_dLat[ix];
- double dLonThis = grid.m_dLon[ix];
-
- double dR =
- sin(dLat0) * sin(dLatThis)
- + cos(dLat0) * cos(dLatThis) * cos(dLonThis - dLon0);
-
- if (dR >= 1.0) {
- dR = 0.0;
- } else if (dR <= -1.0) {
- dR = 180.0;
- } else {
- dR = 180.0 / M_PI * acos(dR);
- }
- if (dR != dR) {
- _EXCEPTIONT("NaN value detected");
- }
-
- if ((ix != ix0) && (dR > dMaxDist)) {
- continue;
- }
-
- // Value at this location
- double dValue = dataState[ix];
-
- // Apply operator
- if (op == ThresholdOp::GreaterThan) {
- if (dValue > dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::LessThan) {
- if (dValue < dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::GreaterThanEqualTo) {
- if (dValue >= dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::LessThanEqualTo) {
- if (dValue <= dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::EqualTo) {
- if (dValue == dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::NotEqualTo) {
- if (dValue != dTargetValue) {
- return true;
- }
-
- } else {
- _EXCEPTIONT("Invalid operation");
- }
-
- // Special case: zero distance
- if (dMaxDist == 0.0) {
- return false;
- }
-
- // Add all neighbors of this point
- for (int n = 0; n < grid.m_vecConnectivity[ix].size(); n++) {
- queueNodes.push(grid.m_vecConnectivity[ix][n]);
- }
- }
-
- return false;
-}
-
-///////////////////////////////////////////////////////////////////////////////
-
// Set of indicator locations stored as grid indices
typedef std::set IndicatorSet;
typedef IndicatorSet::iterator IndicatorSetIterator;
@@ -753,6 +597,12 @@ void ApplyFilters(
_ASSERT(bTag.GetRows() == grid.GetSize());
_ASSERT(bTag.GetRows() == dataState.GetRows());
+ // Threshold operator
+ ThresholdOp threshop;
+ threshop.m_eOp = op;
+ threshop.m_dThresholdValue = dTargetValue;
+ threshop.m_dDistanceDeg = 0.0;
+
// Number of blobs
int nBlobs = 0;
int nBlobsFiltered = 0;
@@ -798,13 +648,8 @@ void ApplyFilters(
// Check if this point satisfies threshold
bool fSatisfiesThreshold =
- SatisfiesThreshold(
- grid,
- dataState,
- iNext,
- op,
- dTargetValue,
- 0.0);
+ threshop.IsSatisfied(
+ grid, dataState, iNext);
if (fSatisfiesThreshold) {
nThresholdPoints++;
@@ -1004,7 +849,6 @@ class DetectBlobsParam {
strTagVar("binary_tag"),
strLongitudeName("lon"),
strLatitudeName("lat"),
- pvecThresholdOp(NULL),
pvecFilterOp(NULL),
pvecGeoFilterOp(NULL),
pvecOutputOp(NULL),
@@ -1057,7 +901,7 @@ class DetectBlobsParam {
std::string strLatitudeName;
// Vector of threshold operators
- std::vector * pvecThresholdOp;
+ ThresholdOpTreeNode aThresholdOp;
// Vector of filter operators
std::vector * pvecFilterOp;
@@ -1088,11 +932,6 @@ void DetectBlobs(
VariableRegistry & varreg,
const DetectBlobsParam & param
) {
- // Dereference pointers to operators
- _ASSERT(param.pvecThresholdOp != NULL);
- std::vector & vecThresholdOp =
- *(param.pvecThresholdOp);
-
// Dereference pointers to operators
_ASSERT(param.pvecFilterOp != NULL);
std::vector & vecFilterOp =
@@ -1265,15 +1104,15 @@ void DetectBlobs(
NcDim * dimTimeOut = NULL;
if ((dimTime != NULL) && (varTime != NULL)) {
if (vecOutputTimes.size() != vecTimes.size()) {
- CopyNcVarTimeSubset(ncInput, ncOutput, "time", vecOutputTimes);
+ CopyNcVarTimeSubset(ncInput, ncOutput, varTime->name(), vecOutputTimes);
} else {
- CopyNcVar(ncInput, ncOutput, "time", true);
+ CopyNcVar(ncInput, ncOutput, varTime->name(), true);
}
dimTimeOut = NcGetTimeDimension(ncOutput);
if (dimTimeOut == NULL) {
- _EXCEPTIONT("Error copying variable \"time\" to output file");
+ _EXCEPTION1("Error copying variable \"%s\" to output file", varTime->name());
}
} /*else if (nAddTimeDim != -1) {
@@ -1456,36 +1295,9 @@ void DetectBlobs(
// Eliminate based on threshold commands
AnnounceStartBlock("Apply threshold commands");
- for (int tc = 0; tc < vecThresholdOp.size(); tc++) {
-
- // Load the search variable data
- Variable & var = varreg.Get(vecThresholdOp[tc].m_varix);
- vecFiles.SetTime(vecTimes[t]);
- var.LoadGridData(varreg, vecFiles, grid);
- const DataArray1D & dataState = var.GetData();
-
- // Loop through data
- for (int i = 0; i < grid.GetSize(); i++) {
- if (bTag[i] == 0) {
- continue;
- }
-
- // Determine if the threshold is satisfied
- bool fSatisfiesThreshold =
- SatisfiesThreshold(
- grid,
- dataState,
- i,
- vecThresholdOp[tc].m_eOp,
- vecThresholdOp[tc].m_dValue,
- vecThresholdOp[tc].m_dDistance
- );
-
- if (!fSatisfiesThreshold) {
- bTag[i] = 0;
- }
- }
- }
+ vecFiles.SetTime(vecTimes[t]);
+ param.aThresholdOp.IsSatisfied(
+ varreg, vecFiles, grid, bTag);
AnnounceEndBlock("Done");
// Eliminate based on filter commands
@@ -1790,30 +1602,9 @@ try {
}
// Parse the threshold operator command string
- std::vector vecThresholdOp;
- dbparam.pvecThresholdOp = &vecThresholdOp;
-
if (strThresholdCmd != "") {
AnnounceStartBlock("Parsing threshold operations");
-
- int iLast = 0;
- for (int i = 0; i <= strThresholdCmd.length(); i++) {
-
- if ((i == strThresholdCmd.length()) ||
- (strThresholdCmd[i] == ';') ||
- (strThresholdCmd[i] == ':')
- ) {
- std::string strSubStr =
- strThresholdCmd.substr(iLast, i - iLast);
-
- int iNextOp = (int)(vecThresholdOp.size());
- vecThresholdOp.resize(iNextOp + 1);
- vecThresholdOp[iNextOp].Parse(varreg, strSubStr);
-
- iLast = i + 1;
- }
- }
-
+ dbparam.aThresholdOp.Parse(varreg, strThresholdCmd);
AnnounceEndBlock("Done");
}
diff --git a/src/nodes/DetectNodes.cpp b/src/nodes/DetectNodes.cpp
index bb56ddb..9ab2d68 100644
--- a/src/nodes/DetectNodes.cpp
+++ b/src/nodes/DetectNodes.cpp
@@ -176,122 +176,6 @@ bool HasClosedContour(
///////////////////////////////////////////////////////////////////////////////
-///
-/// Determine if the given field satisfies the threshold.
-///
-template
-bool SatisfiesThreshold(
- const SimpleGrid & grid,
- const DataArray1D & dataState,
- const int ix0,
- const ThresholdOp::Operation op,
- const double dTargetValue,
- const double dMaxDist
-) {
- // Verify that dMaxDist is less than 180.0
- if (dMaxDist > 180.0) {
- _EXCEPTIONT("MaxDist must be less than 180.0");
- }
-
- // Queue of nodes that remain to be visited
- std::queue queueNodes;
- queueNodes.push(ix0);
-
- // Set of nodes that have already been visited
- std::set setNodesVisited;
-
- // Latitude and longitude at the origin
- double dLat0 = grid.m_dLat[ix0];
- double dLon0 = grid.m_dLon[ix0];
-
- // Loop through all latlon elements
- while (queueNodes.size() != 0) {
- int ix = queueNodes.front();
- queueNodes.pop();
-
- if (setNodesVisited.find(ix) != setNodesVisited.end()) {
- continue;
- }
-
- setNodesVisited.insert(ix);
-
- // Great circle distance to this element
- double dLatThis = grid.m_dLat[ix];
- double dLonThis = grid.m_dLon[ix];
-
- double dR =
- sin(dLat0) * sin(dLatThis)
- + cos(dLat0) * cos(dLatThis) * cos(dLonThis - dLon0);
-
- if (dR >= 1.0) {
- dR = 0.0;
- } else if (dR <= -1.0) {
- dR = 180.0;
- } else {
- dR = 180.0 / M_PI * acos(dR);
- }
- if (dR != dR) {
- _EXCEPTIONT("NaN value detected");
- }
-
- if ((ix != ix0) && (dR > dMaxDist)) {
- continue;
- }
-
- // Value at this location
- double dValue = dataState[ix];
-
- // Apply operator
- if (op == ThresholdOp::GreaterThan) {
- if (dValue > dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::LessThan) {
- if (dValue < dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::GreaterThanEqualTo) {
- if (dValue >= dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::LessThanEqualTo) {
- if (dValue <= dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::EqualTo) {
- if (dValue == dTargetValue) {
- return true;
- }
-
- } else if (op == ThresholdOp::NotEqualTo) {
- if (dValue != dTargetValue) {
- return true;
- }
-
- } else {
- _EXCEPTIONT("Invalid operation");
- }
-
- // Special case: zero distance
- if (dMaxDist == 0.0) {
- return false;
- }
-
- // Add all neighbors of this point
- for (int n = 0; n < grid.m_vecConnectivity[ix].size(); n++) {
- queueNodes.push(grid.m_vecConnectivity[ix][n]);
- }
- }
-
- return false;
-}
-
-///////////////////////////////////////////////////////////////////////////////
-
class DetectCyclonesParam {
public:
@@ -313,7 +197,6 @@ class DetectCyclonesParam {
dMergeDist(0.0),
pvecClosedContourOp(NULL),
pvecNoClosedContourOp(NULL),
- pvecThresholdOp(NULL),
pvecOutputOp(NULL),
nTimeStride(1),
strLatitudeName("lat"),
@@ -367,8 +250,8 @@ class DetectCyclonesParam {
// Vector of no closed contour operators
std::vector * pvecNoClosedContourOp;
- // Vector of threshold operators
- std::vector * pvecThresholdOp;
+ // Threshold operator tree
+ ThresholdOpTreeNode aThresholdOp;
// Vector of output operators
std::vector * pvecOutputOp;
@@ -447,10 +330,6 @@ void DetectCyclonesUnstructured(
std::vector & vecNoClosedContourOp =
*(param.pvecNoClosedContourOp);
- _ASSERT(param.pvecThresholdOp != NULL);
- std::vector & vecThresholdOp =
- *(param.pvecThresholdOp);
-
_ASSERT(param.pvecOutputOp != NULL);
std::vector & vecOutputOp =
*(param.pvecOutputOp);
@@ -617,9 +496,10 @@ void DetectCyclonesUnstructured(
int nRejectedLocation = 0;
int nRejectedMerge = 0;
- DataArray1D vecRejectedClosedContour(vecClosedContourOp.size());
- DataArray1D vecRejectedNoClosedContour(vecNoClosedContourOp.size());
- DataArray1D vecRejectedThreshold(vecThresholdOp.size());
+ std::vector vecRejectedClosedContour(vecClosedContourOp.size(), 0);
+ std::vector vecRejectedNoClosedContour(vecNoClosedContourOp.size(), 0);
+ //std::vector vecRejectedThreshold(param.aThresholdOp.size(), 0);
+ int nRejectedThreshold = 0;
// Eliminate based on interval
if ((param.dMinLatitude != param.dMaxLatitude) ||
@@ -775,41 +655,12 @@ void DetectCyclonesUnstructured(
}
// Eliminate based on thresholds
- for (int tc = 0; tc < vecThresholdOp.size(); tc++) {
-
- std::set setNewCandidates;
-
- // Load the search variable data
- Variable & var = varreg.Get(vecThresholdOp[tc].m_varix);
+ {
vecFiles.SetTime(vecTimes[t]);
- var.LoadGridData(varreg, vecFiles, grid);
- const DataArray1D & dataState = var.GetData();
-
- // Loop through all pressure minima
- std::set::const_iterator iterCandidate
- = setCandidates.begin();
-
- for (; iterCandidate != setCandidates.end(); iterCandidate++) {
-
- // Determine if the threshold is satisfied
- bool fSatisfiesThreshold =
- SatisfiesThreshold(
- grid,
- dataState,
- *iterCandidate,
- vecThresholdOp[tc].m_eOp,
- vecThresholdOp[tc].m_dValue,
- vecThresholdOp[tc].m_dDistance
- );
-
- // If not rejected, add to new pressure minima array
- if (fSatisfiesThreshold) {
- setNewCandidates.insert(*iterCandidate);
- } else {
- vecRejectedThreshold[tc]++;
- }
- }
-
+ std::set setNewCandidates;
+ param.aThresholdOp.IsSatisfied(
+ varreg, vecFiles, grid, setCandidates, setNewCandidates);
+ nRejectedThreshold = static_cast(setCandidates.size() - setNewCandidates.size());
setCandidates = setNewCandidates;
}
@@ -892,20 +743,21 @@ void DetectCyclonesUnstructured(
Announce("Total candidates: %i", setCandidates.size());
Announce("Rejected ( location): %i", nRejectedLocation);
Announce("Rejected ( merged): %i", nRejectedMerge);
-
- for (int tc = 0; tc < vecRejectedThreshold.GetRows(); tc++) {
+ Announce("Rejected ( threshold): %i", nRejectedThreshold);
+/*
+ for (int tc = 0; tc < vecRejectedThreshold.size(); tc++) {
Announce("Rejected (thresh. %s): %i",
- varreg.GetVariableString(vecThresholdOp[tc].m_varix).c_str(),
+ param.aThresholdOp.ToString(varreg, tc).c_str(),
vecRejectedThreshold[tc]);
}
-
- for (int ccc = 0; ccc < vecRejectedClosedContour.GetRows(); ccc++) {
+*/
+ for (int ccc = 0; ccc < vecRejectedClosedContour.size(); ccc++) {
Announce("Rejected (contour %s): %i",
varreg.GetVariableString(vecClosedContourOp[ccc].m_varix).c_str(),
vecRejectedClosedContour[ccc]);
}
- for (int ccc = 0; ccc < vecRejectedNoClosedContour.GetRows(); ccc++) {
+ for (int ccc = 0; ccc < vecRejectedNoClosedContour.size(); ccc++) {
Announce("Rejected (nocontour %s): %i",
varreg.GetVariableString(vecNoClosedContourOp[ccc].m_varix).c_str(),
vecRejectedNoClosedContour[ccc]);
@@ -1230,30 +1082,9 @@ try {
}
// Parse the threshold operator command string
- std::vector vecThresholdOp;
- dcuparam.pvecThresholdOp = &vecThresholdOp;
-
if (strThresholdCmd != "") {
AnnounceStartBlock("Parsing threshold operations");
-
- int iLast = 0;
- for (int i = 0; i <= strThresholdCmd.length(); i++) {
-
- if ((i == strThresholdCmd.length()) ||
- (strThresholdCmd[i] == ';') ||
- (strThresholdCmd[i] == ':')
- ) {
- std::string strSubStr =
- strThresholdCmd.substr(iLast, i - iLast);
-
- int iNextOp = (int)(vecThresholdOp.size());
- vecThresholdOp.resize(iNextOp + 1);
- vecThresholdOp[iNextOp].Parse(varreg, strSubStr);
-
- iLast = i + 1;
- }
- }
-
+ dcuparam.aThresholdOp.Parse(varreg, strThresholdCmd);
AnnounceEndBlock("Done");
}