Commit b0c1efd4 authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Merge branch 'log-interpolation' into 'master'

Log interpolation

See merge request jap/radix!67
parents 60af0a0f 348bc57d
Pipeline #39840 passed with stages
in 15 minutes and 9 seconds
......@@ -60,7 +60,7 @@ linux_gcc_testing:
- module load cmake qt/5.9.0 vtk/8.1.0
- which cmake
- module load valgrind
- cmake -DBUILDNAME=$(uname -s)-GCC-4.8.5-Release-${CI_BUILD_REF_NAME} -DCMAKE_BUILD_TYPE=Release -Dradix_ENABLE_TESTS=ON -Dradix_ENABLE_SECONDARY_TESTED_CODE=ON -Dradix_ENABLE_TESTS=ON -DTPL_ENABLE_VTK=ON -Dradix_ENABLE_radixplot=OFF -Dradix_ENABLE_radixwidgets=OFF ..
- cmake -DBUILDNAME=$(uname -s)-GCC-4.8.5-Release-${CI_BUILD_REF_NAME} -DCMAKE_BUILD_TYPE=Release -Dradix_ENABLE_TESTS=ON -DENABLE_PYTHON_WRAPPERS=ON -Dradix_ENABLE_SECONDARY_TESTED_CODE=ON -Dradix_ENABLE_TESTS=ON -DTPL_ENABLE_VTK=ON -Dradix_ENABLE_radixplot=OFF -Dradix_ENABLE_radixwidgets=OFF ..
- ctest -D ExperimentalStart -D ExperimentalBuild -D ExperimentalTest -D ExperimentalSubmit
linux_analysis:
......
......@@ -419,29 +419,29 @@ int main(int argc, char **argv)
}
std::cout << " Temperature:" << std::endl;
inputTemps = interpolateToOtherBaseValues(inputPressures,
resampledPressures, inputTemps);
inputTemps = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputTemps, true);
if (usingRelHum)
{
std::cout << " Relative humidity:" << std::endl;
inputRelHums = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputRelHums);
inputPressures, resampledPressures, inputRelHums, true);
}
if (usingDewPt)
{
std::cout << " Dew point:" << std::endl;
inputDewPts = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputDewPts);
inputPressures, resampledPressures, inputDewPts, true);
}
std::cout << " Wind speed:" << std::endl;
inputWSpds = interpolateToOtherBaseValues(inputPressures,
resampledPressures, inputWSpds);
inputWSpds = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputWSpds, true);
std::cout << " Wind direction:" << std::endl;
inputWDirs = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputWDirs, true);
inputPressures, resampledPressures, inputWDirs, true, true);
std::cout << " Height:" << std::endl;
inputHeights = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputHeights);
inputPressures, resampledPressures, inputHeights, true);
std::cout << " Replacing pressure data with resampled array" << std::endl;
radix_block(
......@@ -457,23 +457,23 @@ int main(int argc, char **argv)
{
std::cout << " Not resampling data" << std::endl;
std::cout << " Temperature:" << std::endl;
inputTemps = interpolateValues(inputPressures, inputTemps);
inputTemps = interpolateValues(inputPressures, inputTemps, true);
if (usingRelHum)
{
std::cout << " Relative humidity:" << std::endl;
inputRelHums = interpolateValues(inputPressures, inputRelHums);
inputRelHums = interpolateValues(inputPressures, inputRelHums, true);
}
if (usingDewPt)
{
std::cout << " Dew point:" << std::endl;
inputDewPts = interpolateValues(inputPressures, inputDewPts);
inputDewPts = interpolateValues(inputPressures, inputDewPts, true);
}
std::cout << " Wind speed:" << std::endl;
inputWSpds = interpolateValues(inputPressures, inputWSpds);
inputWSpds = interpolateValues(inputPressures, inputWSpds, false);
std::cout << " Wind direction:" << std::endl;
inputWDirs = interpolateValues(inputPressures, inputWDirs, true);
inputWDirs = interpolateValues(inputPressures, inputWDirs, false, true);
std::cout << " Height:" << std::endl;
inputHeights = interpolateValues(inputPressures, inputHeights);
inputHeights = interpolateValues(inputPressures, inputHeights, true);
}
std::cout << " Interpolation complete." << std::endl;
......
......@@ -15,13 +15,14 @@ namespace radix
template <typename T>
std::vector<T> RADIX_PUBLIC interpolateValues(
const std::vector<T> &baseValues, const std::vector<T> &valuesToInterpolate,
const bool circular = false, const T missingValue = T(-9999));
const bool logarithmic = false, const bool circular = false,
const T missingValue = T(-9999));
template <typename T>
std::vector<T> RADIX_PUBLIC interpolateToOtherBaseValues(
const std::vector<T> &baseValues, const std::vector<T> &newBaseValues,
const std::vector<T> &valuesToInterpolate, const bool circular = false,
const T missingValues = T(-9999));
const std::vector<T> &valuesToInterpolate, const bool logarithmic = false,
const bool circular = false, const T missingValue = T(-9999));
} // namespace radix
......
......@@ -15,6 +15,8 @@ namespace radix
* to be present. Used for determining 'interpolation bounds' (there's
* definitely a better phrase for this)
* @param valuesToInterpolate The values to interpolate
* @param logarithmic Flag indicating whether to interpolate linearly with
* baseValues (false) or by their natural log (true)
* @param circular Flag indicating whether to use circular interpolation (for
* angles in degrees)
* @param missingValue The value which indicates that the value is missing
......@@ -22,7 +24,8 @@ namespace radix
template <typename T>
std::vector<T> interpolateValues(const std::vector<T> &baseValues,
const std::vector<T> &valuesToInterpolate,
const bool circular, const T missingValue)
const bool logarithmic, const bool circular,
const T missingValue)
{
std::vector<T> interpolatedValues = valuesToInterpolate;
......@@ -109,10 +112,21 @@ std::vector<T> interpolateValues(const std::vector<T> &baseValues,
}
}
interpolatedValues[j] =
lastGoodValue + ((nextGoodValue - lastGoodValue) *
((baseValues[j] - baseValues[lastGood]) /
(baseValues[nextGood] - baseValues[lastGood])));
if (logarithmic)
{
interpolatedValues[j] =
lastGoodValue +
((nextGoodValue - lastGoodValue) *
((log(baseValues[j]) - log(baseValues[lastGood])) /
(log(baseValues[nextGood]) - log(baseValues[lastGood]))));
}
else
{
interpolatedValues[j] =
lastGoodValue + ((nextGoodValue - lastGoodValue) *
((baseValues[j] - baseValues[lastGood]) /
(baseValues[nextGood] - baseValues[lastGood])));
}
if (circular)
{
......@@ -128,18 +142,22 @@ std::vector<T> interpolateValues(const std::vector<T> &baseValues,
/**
* @brief interpolateToOtherBaseValues
* @param baseValues
* @param newBaseValues
* @param valuesToInterpolate
* @param circular
* @param missingValues
* @param baseValues Base vector for interpolation - requires all values
* to be present. Used for determining 'interpolation bounds'
* @param newBaseValues The new vector of base values to use
* @param valuesToInterpolate The values to interpolate to the new base values
* @param logarithmic Flag indicating whether to interpolate linearly with
* baseValues (false) or by their natural log (true)
* @param circular Flag indicating whether to use circular interpolation (for
* angles in degrees)
* @param missingValue The value which indicates that the value is missing
* @return
*/
template <typename T>
std::vector<T> interpolateToOtherBaseValues(
const std::vector<T> &baseValues, const std::vector<T> &newBaseValues,
const std::vector<T> &valuesToInterpolate, const bool circular,
const T missingValues)
const std::vector<T> &valuesToInterpolate, const bool logarithmic,
const bool circular, const T missingValue)
{
std::vector<T> finalInterpolatedValues(newBaseValues.size());
......@@ -148,7 +166,7 @@ std::vector<T> interpolateToOtherBaseValues(
// First interpolate across the existing data to ensure there are no missing
radix_line(" Initial interpolation to remove missing values in input data");
std::vector<T> initialInterpolatedValues = interpolateValues(
baseValues, valuesToInterpolate, circular, missingValues);
baseValues, valuesToInterpolate, logarithmic, circular, missingValue);
// If this fails, we have to fail this one too
if (initialInterpolatedValues.size() == 0)
{
......@@ -250,11 +268,23 @@ std::vector<T> interpolateToOtherBaseValues(
}
else
{
finalInterpolatedValues[i] =
lastGoodValue +
((nextGoodValue - lastGoodValue) *
((newBaseValues[i] - baseValues[lowerBaseIndex]) /
(baseValues[higherBaseIndex] - baseValues[lowerBaseIndex])));
if (logarithmic)
{
finalInterpolatedValues[i] =
lastGoodValue +
((nextGoodValue - lastGoodValue) *
((log(newBaseValues[i]) - log(baseValues[lowerBaseIndex])) /
(log(baseValues[higherBaseIndex]) -
log(baseValues[lowerBaseIndex]))));
}
else
{
finalInterpolatedValues[i] =
lastGoodValue +
((nextGoodValue - lastGoodValue) *
((newBaseValues[i] - baseValues[lowerBaseIndex]) /
(baseValues[higherBaseIndex] - baseValues[lowerBaseIndex])));
}
}
if (circular)
......
MESSAGE(STATUS "${PACKAGE_NAME} enabling swig python bindings.")
INCLUDE(RadixSWIGPyModules)
......
......@@ -14,10 +14,10 @@ namespace radix
{
template<typename T>
std::vector<T> interpolateValues(std::vector<T> const & baseValues, std::vector<T> const & valuesToInterpolate, const bool circular, const T missingValue);
std::vector<T> interpolateValues(std::vector<T> const & baseValues, std::vector<T> const & valuesToInterpolate, const bool logarithmic, const bool circular, const T missingValue);
template<typename T>
std::vector<T> interpolateToOtherBaseValues(std::vector<T> const & baseValues, std::vector<T> const & newBaseValues, std::vector<T> const & valuesToInterpolate, const bool circular, const T missingValue);
std::vector<T> interpolateToOtherBaseValues(std::vector<T> const & baseValues, std::vector<T> const & newBaseValues, std::vector<T> const & valuesToInterpolate, const bool logarithmic, const bool circular, const T missingValue);
%template(interpolate_values) interpolateValues<double>;
%template(interpolate_to_other_base_values) interpolateToOtherBaseValues<double>;
......
......@@ -4,3 +4,18 @@ ADD_GOOGLE_TEST(tstInterpolate.cc NP 1)
ADD_GOOGLE_TEST(tstLognormal.cc NP 1)
ADD_GOOGLE_TEST(tstMatrix.cc NP 1)
ADD_GOOGLE_TEST(tstUtil.cc NP 1)
# Add the Python test
IF(ENABLE_PYTHON_WRAPPERS)
SET(INTERPOLATE_PYTHON_TEST_NAME tst_interpolate_py)
SET(INTERPOLATE_PYTHON_TEST_FILE test_interpolate.py)
TRIBITS_ADD_ADVANCED_TEST(
${INTERPOLATE_PYTHON_TEST_NAME}
OVERALL_WORKING_DIRECTORY TEST_NAME
TEST_0
CMND ${PYTHON_EXECUTABLE}
ARGS "${CMAKE_CURRENT_SOURCE_DIR}/${INTERPOLATE_PYTHON_TEST_FILE}"
PASS_REGULAR_EXPRESSION "OK"
ENVIRONMENT "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}/../python/"
)
ENDIF()
......@@ -10,30 +10,76 @@ circular = False
class TestInterpolateBindings(unittest.TestCase):
def test_interpolate_basic(self):
# Basic test of interpolate_values
print("Testing interpolate_values...")
print("Testing interpolate_values (linear)...")
logarithmic = False
base_values = (10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0)
interp_values = (1.0, missing_value, 3.0, missing_value, 6.0,
missing_value, 1.0)
expect_values = (1.0, 2.0, 3.0, 4.5, 6.0, 3.5, 1.0)
test_values = radixmath.interpolate_values(
base_values, interp_values, circular, missing_value
base_values, interp_values, logarithmic, circular, missing_value
)
for i in range(0, len(expect_values)):
self.assertAlmostEqual(expect_values[i], test_values[i])
def test_interpolate_log_basic(self):
# Basic test of interpolate_values
print("Testing interpolate_values (logarithmic)...")
logarithmic = True
base_values = (1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0)
interp_values = (1.0, missing_value, 3.0, missing_value, 6.0,
missing_value, 1.0)
expect_values = (1.0, 2.0, 3.0, 4.5, 6.0, 3.5, 1.0)
test_values = radixmath.interpolate_values(
base_values, interp_values, logarithmic, circular, missing_value
)
self.assertEqual(expect_values, test_values)
for i in range(0, len(expect_values)):
self.assertAlmostEqual(expect_values[i], test_values[i])
def test_interpolate_to_other_base_basic(self):
print("Testing interpolate_to_other_base_values")
print("Testing interpolate_to_other_base_values (linear)...")
logarithmic = False
base_values = (1000.0, 2000.0, 3000.0)
new_base_values = (500.0, 1500.0, 2500.0, 3500.0)
interp_values = (10.5, missing_value, 30.5)
expect_values = (10.5, 15.5, 25.5, 30.5)
test_values = radixmath.interpolate_to_other_base_values(
base_values, new_base_values, interp_values, circular, missing_value
base_values, new_base_values, interp_values, logarithmic,
circular, missing_value
)
for i in range(0, len(expect_values)):
self.assertAlmostEqual(expect_values[i], test_values[i])
def test_interpolate_to_other_base_basic_log(self):
print("Testing interpolate_to_other_base_values (logarithmic)...")
logarithmic = True
base_values = (1000.0, 4000.0, 16000.0)
new_base_values = (500.0, 2000.0, 8000.0, 32000.0)
interp_values = (10.5, missing_value, 30.5)
expect_values = (10.5, 15.5, 25.5, 30.5)
test_values = radixmath.interpolate_to_other_base_values(
base_values, new_base_values, interp_values, logarithmic,
circular, missing_value
)
self.assertEqual(expect_values, test_values)
for i in range(0, len(expect_values)):
self.assertAlmostEqual(expect_values[i], test_values[i])
if __name__ == '__main__':
unittest.main()
......@@ -6,7 +6,7 @@
using namespace radix;
TEST(Interpolate, InternalInterpolation)
TEST(Interpolate, InternalInterpolationLinear)
{
float missingValue = -9999.f;
float tolerance = 0.0001f;
......@@ -42,7 +42,38 @@ TEST(Interpolate, InternalInterpolation)
std::vector<float> interpValues3{16.f, missingValue, missingValue, 2.f};
std::vector<float> expectValues3{16.f, 8.f, 4.f, 2.f};
std::vector<float> testValues3 =
interpolateValues(baseValues3, interpValues3, false, missingValue);
interpolateValues(baseValues3, interpValues3, false, false, missingValue);
ASSERT_EQ(expectValues3.size(), testValues3.size());
for (size_t i = 0; i < testValues3.size(); ++i)
{
EXPECT_NEAR(expectValues3[i], testValues3[i], expectValues3[i] * tolerance);
}
}
TEST(Interpolate, InternalInterpolationLog)
{
float missingValue = -9999.f;
float tolerance = 0.0001f;
std::vector<float> baseValues1{1.f, 10.f, 100.f, 1000.f,
10000.f, 100000.f, 1000000.f};
std::vector<float> interpValues1{1.f, missingValue, 3.f, missingValue,
6.f, missingValue, 1.f};
std::vector<float> expectValues1{1.f, 2.f, 3.f, 4.5, 6.f, 3.5, 1.f};
std::vector<float> testValues1 =
interpolateValues(baseValues1, interpValues1, true);
ASSERT_EQ(expectValues1.size(), testValues1.size());
for (size_t i = 0; i < testValues1.size(); ++i)
{
EXPECT_NEAR(expectValues1[i], testValues1[i], expectValues1[i] * tolerance);
}
missingValue = 3e12f;
std::vector<float> baseValues3{800.f, 400.f, 200.f, 100.f};
std::vector<float> interpValues3{8.f, missingValue, missingValue, 2.f};
std::vector<float> expectValues3{8.f, 6.f, 4.f, 2.f};
std::vector<float> testValues3 =
interpolateValues(baseValues3, interpValues3, true, false, missingValue);
ASSERT_EQ(expectValues3.size(), testValues3.size());
for (size_t i = 0; i < testValues3.size(); ++i)
{
......@@ -81,7 +112,7 @@ TEST(Interpolate, CircularInterpolation)
350.f, missingValue, 30.f};
std::vector<float> expectValues1{90.f, 180.f, 270.f, 350.f, 10.f, 30.f};
std::vector<float> testValues1 =
interpolateValues(baseValues1, interpValues1, true, missingValue);
interpolateValues(baseValues1, interpValues1, false, true, missingValue);
ASSERT_EQ(expectValues1.size(), testValues1.size());
for (size_t i = 0; i < testValues1.size(); ++i)
{
......@@ -164,6 +195,47 @@ TEST(Interpolate, NewBaseValues)
}
}
TEST(Interpolate, NewBaseValuesLog)
{
float missingValue = -9999.f;
float tolerance = 0.0001f;
std::vector<float> baseValues1{1000.f, 4000.f, 16000.f};
std::vector<float> newBaseValues1{500.f, 2000.f, 8000.f, 32000.f};
std::vector<float> interpValues1{10.5f, 20.5f, 30.5f};
std::vector<float> expectValues1{10.5f, 15.5f, 25.5f, 30.5f};
std::vector<float> testValues1 = interpolateToOtherBaseValues(
baseValues1, newBaseValues1, interpValues1, true);
ASSERT_EQ(expectValues1.size(), testValues1.size());
for (size_t i = 0; i < testValues1.size(); ++i)
{
EXPECT_NEAR(expectValues1[i], testValues1[i], expectValues1[i] * tolerance);
}
std::vector<float> interpValuesMissing{10.5f, missingValue, 30.5f};
std::vector<float> testValuesMissing = interpolateToOtherBaseValues(
baseValues1, newBaseValues1, interpValuesMissing, true);
ASSERT_EQ(expectValues1.size(), testValuesMissing.size());
for (size_t i = 0; i < testValuesMissing.size(); ++i)
{
EXPECT_NEAR(expectValues1[i], testValuesMissing[i],
expectValues1[i] * tolerance);
}
// High-to-low base values
std::vector<float> baseValues2{16000.f, 4000.f, 1000.f};
std::vector<float> newBaseValues2{32000.f, 8000.f, 2000.f, 500.f};
std::vector<float> interpValues2{30.5f, 20.5f, 10.5f};
std::vector<float> expectValues2{30.5f, 25.5f, 15.5f, 10.5f};
std::vector<float> testValues2 = interpolateToOtherBaseValues(
baseValues2, newBaseValues2, interpValues2, true);
ASSERT_EQ(expectValues2.size(), testValues2.size());
for (size_t i = 0; i < testValues2.size(); ++i)
{
EXPECT_NEAR(expectValues2[i], testValues2[i], expectValues2[i] * tolerance);
}
}
TEST(Interpolate, NewBaseValuesDoubles)
{
double missingValue = -9999.0;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment