Temporary fix for missing values/NaN incompatibility for SyCLoPS

README.md

@@ -27,7 +27,7 @@ To install from conda use the command line:
 
 # Installation via CMake (Recommended)
 
-TempestExtremes can be built and installed on various systems using CMake. Our new script, `./quick_make_unix.sh`, automatically detects your platform(UNIX-based systems only) and loads any required modules before building.
+TempestExtremes can be built and installed on various systems using CMake. Our new script, `./quick_make_unix.sh` in the TempestExtremes root directory, automatically detects your operating system (Unix-based systems only) and loads any required modules before building. To install all TempestExtremes executables to the ./bin directory, just run the quick make script.
 
 ## General CMake Configuration
 

src/base/DataOp.cpp

@@ -436,9 +436,14 @@ bool DataOp_VECMAG::Apply(
 	const DataArray1D<float> & dataRight = *(vecArgData[1]);
 
 	for (int i = 0; i < dataout.GetRows(); i++) {
-		dataout[i] =
-			sqrt(dataLeft[i] * dataLeft[i]
-				+ dataRight[i] * dataRight[i]);
+		bool fIsFillValue = dataLeft.IsFillValueAtIx(i);
+		if (!fIsFillValue){
+			dataout[i] =
+				sqrt(dataLeft[i] * dataLeft[i]
+					+ dataRight[i] * dataRight[i]);
+		} else {
+				dataout[i]=std::numeric_limits<double>::quiet_NaN();
+				}
 	}
 
 	return true;
@@ -699,7 +704,12 @@ bool DataOp_DIFF::Apply(
 		const DataArray1D<float> & dataRight = *(vecArgData[1]);
 
 		for (int i = 0; i < dataout.GetRows(); i++) {
-			dataout[i] = dataLeft[i] - dataRight[i];
+			bool fIsFillValue = (dataLeft.IsFillValueAtIx(i) || dataRight.IsFillValueAtIx(i));
+			if (!fIsFillValue){
+				dataout[i] = dataLeft[i] - dataRight[i];
+			} else {
+				dataout[i]=std::numeric_limits<double>::quiet_NaN();
+			}
 		}
 	}
 
@@ -1702,7 +1712,6 @@ void BuildCurlOperator(
 ) {
 	opCurlE.Clear();
 	opCurlN.Clear();
-
 	int iRef = 0;
 
 	// Scaling factor used in Curl calculation
@@ -1955,7 +1964,6 @@ bool DataOp_CURL::Apply(
 
 		m_fInitialized = true;
 	}
-
 	m_opCurlE.Apply(*(vecArgData[0]), dataout, true);
 	m_opCurlN.Apply(*(vecArgData[1]), dataout, false); 
 

src/base/SparseMatrix.h

@@ -18,7 +18,6 @@
 #define _SPARSEMATRIX_H_
 
 #include "DataArray1D.h"
-
 #include <map>
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -190,8 +189,13 @@ class SparseMatrix {
 
 		SparseMapConstIterator iter = m_mapEntries.begin();
 		for (; iter != m_mapEntries.end(); iter++) {
-			dataVectorOut[iter->first.first] +=
+			if (dataVectorIn[iter->first.second]>1e14){
+				// dataVectorOut[iter->first.first] += std::numeric_limits<double>::quiet_NaN();
+				continue;
+			} else {
+				dataVectorOut[iter->first.first] +=
 				iter->second * dataVectorIn[iter->first.second];
+			}
 		}
 	}
 

src/nodes/NodeFileFilter.cpp

@@ -69,7 +69,7 @@ class NearbyBlobsOp {
 		m_varix(InvalidVariableIndex),
 		m_dDistanceDeg(0.0),
 		m_eOp(GreaterThan),
-		m_dValue(0.0),
+		m_strValue(""),
 		m_dMaxDistDeg(180.0)
 	{ }
 
@@ -133,7 +133,7 @@ class NearbyBlobsOp {
 
 				// Read in value
 				} else if (eReadMode == ReadMode_Value) {
-					m_dValue = atof(strSubStr.c_str());
+					m_strValue = strSubStr;
 
 					iLast = i + 1;
 					eReadMode = ReadMode_MaxDist;
@@ -191,12 +191,7 @@ class NearbyBlobsOp {
 			strDescription += " is not equal to ";
 		}
 
-		if (fabs(m_dValue) < 1.0e-4) {
-			snprintf(szBuffer, 128, "%e", m_dValue);
-		} else {
-			snprintf(szBuffer, 128, "%f", m_dValue);
-		}
-		strDescription += szBuffer;
+		strDescription += m_strValue;
 
 		snprintf(szBuffer, 128, "%f", m_dMaxDistDeg);
 		strDescription += std::string(" (max dist ") + szBuffer + " degrees)";
@@ -209,35 +204,36 @@ class NearbyBlobsOp {
 	///		Determine if a particular value satisfies this threshold.
 	///	</summary>
 	bool SatisfiedBy(
-		double dValue
+		double dValue,
+        double dOpValue
 	) const {
 		if (m_eOp == GreaterThan) {
-			if (dValue > m_dValue) {
+			if (dValue > dOpValue) {
 				return true;
 			}
 
 		} else if (m_eOp == LessThan) {
-			if (dValue < m_dValue) {
+			if (dValue < dOpValue) {
 				return true;
 			}
 
 		} else if (m_eOp == GreaterThanEqualTo) {
-			if (dValue >= m_dValue) {
+			if (dValue >= dOpValue) {
 				return true;
 			}
 
 		} else if (m_eOp == LessThanEqualTo) {
-			if (dValue <= m_dValue) {
+			if (dValue <= dOpValue) {
 				return true;
 			}
 
 		} else if (m_eOp == EqualTo) {
-			if (dValue == m_dValue) {
+			if (dValue == dOpValue) {
 				return true;	
 			}
 
 		} else if (m_eOp == NotEqualTo) {
-			if (dValue != m_dValue) {
+			if (dValue != dOpValue) {
 				return true;
 			}
 
@@ -267,7 +263,7 @@ class NearbyBlobsOp {
 	///	<summary>
 	///		Threshold value.
 	///	</summary>
-	double m_dValue;
+	std::string m_strValue;
 
 	///	<summary>
 	///		Maximum distance in degrees.
@@ -441,6 +437,7 @@ void BuildMask_ByContour(
 template <typename real>
 void BuildMask_NearbyBlobs(
 	const SimpleGrid & grid,
+    const NodeFile & nodefile,
 	const DataArray1D<real> & dataState,
 	std::vector<double> & vecLonRad,
 	std::vector<double> & vecLatRad,
@@ -450,6 +447,16 @@ void BuildMask_NearbyBlobs(
 	// Get the variable
 	const double dDistDeg = nearbyblobsop.m_dDistanceDeg;
 	const double dMaxDistDeg = nearbyblobsop.m_dMaxDistDeg;
+    std::string strValue = nearbyblobsop.m_strValue;
+    double dOpValue = 0.0;
+
+    std::vector<double> VecValue;
+    if (STLStringHelper::IsFloat(strValue)) {
+        dOpValue = std::stod(strValue);
+    } else {
+        nodefile.InterpolatedColumnDouble(
+			strValue, VecValue);
+    }
 
 	_ASSERT(vecLonRad.size() == vecLatRad.size());
 	_ASSERT(dataState.GetRows() == grid.GetSize());
@@ -472,7 +479,7 @@ void BuildMask_NearbyBlobs(
 		std::queue<int> queueNodes;
 		queueNodes.push(ix0);
 
-		// Set of nodes that have already beenvisited
+		// Set of nodes that have already been visited
 		std::set<int> setNodesVisited;
 
 		// Loop through all elements
@@ -504,8 +511,12 @@ void BuildMask_NearbyBlobs(
 				queueNodes.push(grid.m_vecConnectivity[ix][n]);
 			}
 
+            if (!VecValue.empty()){
+                dOpValue=VecValue[j];
+            }
+
 			// Check if this point satisfies the nearbyblobs criteria
-			if (!nearbyblobsop.SatisfiedBy(static_cast<double>(dataState[ix]))) {
+			if (!nearbyblobsop.SatisfiedBy(static_cast<double>(dataState[ix]),dOpValue)) {
 				continue;
 			}
 
@@ -543,7 +554,7 @@ void BuildMask_NearbyBlobs(
 				}
 
 				// Verify this point satisfies the condition
-				if (!nearbyblobsop.SatisfiedBy(static_cast<double>(dataState[ixblob]))) {
+				if (!nearbyblobsop.SatisfiedBy(static_cast<double>(dataState[ixblob]),dOpValue)) {
 
 					// Isn't part of the blob, but add it to the list of
 					// nodes to visit.
@@ -562,9 +573,9 @@ void BuildMask_NearbyBlobs(
 					queueThresholdedNodes.push(grid.m_vecConnectivity[ixblob][n]);
 				}
 
-				dataMask[ixblob] = 1.0;
+				dataMask[ixblob] += 1.0;
 			}
-		}
+		}			
 	}
 }
 
@@ -976,23 +987,36 @@ void NodeFileFilter(
 
 			AnnounceEndBlock("Done");
 		}
+		int tc = 0;
+
 		if (vecNearbyBlobsOp.size() != 0) {
 			AnnounceStartBlock("Building mask (nearbyblobs)");
-			Variable & varOp = varreg.Get(vecNearbyBlobsOp[0].m_varix);
-			vecInFiles.SetConstantTimeIx(t);
-			varOp.LoadGridData(varreg, vecInFiles, grid);
-			const DataArray1D<float> & dataState = varOp.GetData();
-			_ASSERT(dataState.GetRows() == grid.GetSize());
-
-			BuildMask_NearbyBlobs<float>(
-				grid,
-				dataState,
-				vecLonRad,
-				vecLatRad,
-				vecNearbyBlobsOp[0],
-				dataMask);
-
-			AnnounceEndBlock("Done");
+			for (tc = 0; tc < vecNearbyBlobsOp.size(); tc++) {
+				Variable & varOp = varreg.Get(vecNearbyBlobsOp[tc].m_varix);
+				vecInFiles.SetConstantTimeIx(t);
+				varOp.LoadGridData(varreg, vecInFiles, grid);
+				const DataArray1D<float> & dataState = varOp.GetData();
+				_ASSERT(dataState.GetRows() == grid.GetSize());
+				int vecsize = vecNearbyBlobsOp.size();
+				BuildMask_NearbyBlobs<float>(
+					grid,
+					nodefile,
+					dataState,
+					vecLonRad,
+					vecLatRad,
+					vecNearbyBlobsOp[tc],
+					dataMask);
+
+				AnnounceEndBlock("Done");
+			}
+        }
+
+		if (tc > 1){
+			for (int i = 0; i < dataMask.GetRows(); i++) {
+				if(dataMask[i] > 0){
+					dataMask[i] = dataMask[i] - (tc-1);
+				}
+			}
 		}
 
 		if (fInvert) {
@@ -1381,9 +1405,23 @@ try {
 
 	if (strNearbyBlobs != "") {
 		AnnounceStartBlock("Parsing --nearbyblobs");
-		NearbyBlobsOp op;
-		op.Parse(varreg, strNearbyBlobs);
-		vecNearbyBlobsOp.push_back(op);
+        int iLast = 0;
+		for (int i = 0; i <= strNearbyBlobs.length(); i++) {
+
+			if ((i == strNearbyBlobs.length()) ||
+				(strNearbyBlobs[i] == ';') ||
+				(strNearbyBlobs[i] == ':')
+			) {
+				std::string strSubStr =
+					strNearbyBlobs.substr(iLast, i - iLast);
+
+				int iNextOp = (int)(vecNearbyBlobsOp.size());
+				vecNearbyBlobsOp.resize(iNextOp + 1);
+				vecNearbyBlobsOp[iNextOp].Parse(varreg, strSubStr);
+
+				iLast = i + 1;
+			}
+		}
 		AnnounceEndBlock(NULL);
 	}