diff --git a/include/schaapcommon/h5parm/soltab.h b/include/schaapcommon/h5parm/soltab.h
index 4860ae623f4d190ce492edff681522913222d7fa..bf93ed6ed29bea007f07758a06e812a1f0428620 100644
--- a/include/schaapcommon/h5parm/soltab.h
+++ b/include/schaapcommon/h5parm/soltab.h
@@ -29,6 +29,30 @@ struct TimesAndFrequencies {
size_t num_freqs;
};
+/**
+ * Determine the axis values (e.g. times) that fall within a requested
+ * range. The unit tests give a number of examples which may help
+ * understanding.
+ * @param axis a sorted non-empty vector with axis values. This vector is
+ * reused as return value to avoid an allocation.
+ * @param requested_start First axis value that is needed.
+ * @param requested_end Last axis value that is needed.
+ * @param nearest If @c true, the requested values are rounded to the nearest
+ * axis values. If @c false, it is assumed that interpolation is required.
+ * This makes sure that there are samples available at the beginning and the
+ * end of the result for interpolation. Particularly, the first selected
+ * value is equal or lower than @p requested_start, and the last
+ * selected value is equal or higher than @p requested_end, if such values
+ * exist. Otherwise, the most extreme values are returned.
+ * @param [out] start_index index of the first axis value returned.
+ * @returns sorted vector with axis values that correspond with the requested
+ * selection.
+ */
+std::vector<double> SelectAxisValues(std::vector<double>&& axis,
+ double requested_start,
+ double requested_end, bool nearest,
+ size_t& start_index);
+
/// @brief SolTab is a solution table as defined in the H5Parm standard. It
/// contains one solution, e.g. all TEC values, with different axes
/// for that solution (e.g. time, freq, pol).
diff --git a/src/h5parm/jonesparameters.cc b/src/h5parm/jonesparameters.cc
index 438f94f0aca2327048c562ec289d52fba2cfd81b..119c4c8a5077bd21f961771c97069eb31b202fff 100644
--- a/src/h5parm/jonesparameters.cc
+++ b/src/h5parm/jonesparameters.cc
@@ -79,6 +79,7 @@ JonesParameters::JonesParameters(
if (parm_size == 0U) {
parm_size = GetNParmValues(gain_type);
}
+ // Indexed as parm_values[parameter][antenna][time-frequency index]
std::vector<std::vector<std::vector<double>>> parm_values(parm_size);
for (auto& parm_values : parm_values) {
parm_values.resize(antenna_names.size());
diff --git a/src/h5parm/soltab.cc b/src/h5parm/soltab.cc
index 0826c840518a91c82d34749a195c3635f7da3e8b..36e568316b88ba2b85599585c78deb8ed2bc8eca 100644
--- a/src/h5parm/soltab.cc
+++ b/src/h5parm/soltab.cc
@@ -39,6 +39,51 @@ std::vector<std::string> Tokenize(const std::string& str,
}
} // namespace
+std::vector<double> SelectAxisValues(std::vector<double>&& axis,
+ double requested_start,
+ double requested_end, bool nearest,
+ size_t& start_index) {
+ assert(!axis.empty());
+ assert(std::is_sorted(axis.cbegin(), axis.cend()));
+ assert(requested_start <= requested_end);
+ // ReadAxes() already checked that the time axis is sorted.
+ std::vector<double>::const_iterator lower =
+ std::lower_bound(axis.cbegin(), axis.cend(), requested_start);
+ std::vector<double>::const_iterator upper;
+ if (nearest) {
+ upper = std::upper_bound(lower, axis.cend(), requested_end);
+ // move upper to nearest element
+ if (upper > axis.cbegin() && upper != axis.cend()) {
+ if (requested_end - *(upper - 1) < *upper - requested_end) {
+ --upper;
+ }
+ }
+ // move lower to the nearest element
+ if (lower > axis.cbegin()) {
+ if (lower == axis.cend() ||
+ requested_start - *(lower - 1) < *lower - requested_start) {
+ --lower;
+ }
+ }
+ } else {
+ upper = std::lower_bound(lower, axis.cend(), requested_end);
+ // lower may be larger than the requested value; if this case move lower
+ // one value back.
+ if (lower > axis.cbegin() &&
+ (lower == axis.cend() || *lower > requested_start)) {
+ --lower;
+ }
+ }
+ // upper must point to the element after the element that is nearest/bounded:
+ if (upper != axis.cend()) {
+ ++upper;
+ }
+ start_index = lower - axis.cbegin();
+ axis.assign(lower, upper);
+ assert(!axis.empty());
+ return std::move(axis);
+}
+
SolTab::SolTab(H5::Group group, const std::string& type,
const std::vector<AxisInfo>& axes)
: H5::Group(group), type_(type), axes_(axes) {
@@ -333,64 +378,25 @@ TimesAndFrequencies SolTab::GetTimesAndFrequencies(
assert(!times.empty());
assert(!frequencies.empty());
- TimesAndFrequencies times_and_frequencies;
- times_and_frequencies.time_axis = {0.0};
- times_and_frequencies.start_time_index = 0;
- times_and_frequencies.num_times = 1;
- times_and_frequencies.freq_axis = {0.0};
- times_and_frequencies.start_freq_index = 0;
- times_and_frequencies.num_freqs = 1;
+ TimesAndFrequencies tf;
+ tf.time_axis = {0.0};
+ tf.start_time_index = 0;
+ tf.freq_axis = {0.0};
+ tf.start_freq_index = 0;
+ tf.num_freqs = 1;
if (HasAxis("time")) {
- times_and_frequencies.time_axis = GetRealAxis("time");
- if (times.size() > 1) {
- times_and_frequencies.num_times = times_and_frequencies.time_axis.size();
- } else {
- // When only one time is requested, find the nearest H5 time indices
- // and only read H5 data for those indices.
-
- // ReadAxes() already checked that the "time" axis is sorted.
- const auto lower = std::lower_bound(
- times_and_frequencies.time_axis.begin(),
- times_and_frequencies.time_axis.end(), times.front());
- if (lower == times_and_frequencies.time_axis.begin()) {
- // times_and_frequencies.start_time_index remains 0
- times_and_frequencies.time_axis.resize(1); // Keep the first item only.
- } else if (lower == times_and_frequencies.time_axis.end()) {
- times_and_frequencies.start_time_index =
- times_and_frequencies.time_axis.size() - 1;
- times_and_frequencies.time_axis = {
- times_and_frequencies.time_axis.back()};
- } else {
- // The time lies between *(lower-1) and *lower.
- times_and_frequencies.start_time_index =
- std::distance(times_and_frequencies.time_axis.begin(), lower) - 1;
- if (nearest) {
- // Only use the nearest entry.
- if (times.front() - *(lower - 1) > *lower - times.front()) {
- ++times_and_frequencies.start_time_index;
- }
- times_and_frequencies.time_axis = {
- times_and_frequencies
- .time_axis[times_and_frequencies.start_time_index]};
- } else {
- // Only load the two entries around the entry.
- times_and_frequencies.num_times = 2;
- times_and_frequencies.time_axis = {*(lower - 1), *lower};
- }
- }
- }
+ tf.time_axis = SelectAxisValues(GetRealAxis("time"), times.front(),
+ times.back(), nearest, tf.start_time_index);
}
+ tf.num_times = tf.time_axis.size();
if (HasAxis("freq")) {
std::vector<double> full_freq_axis_h5 = GetRealAxis("freq");
- times_and_frequencies.start_freq_index = GetFreqIndex(frequencies.front());
- times_and_frequencies.num_freqs = GetFreqIndex(frequencies.back()) -
- times_and_frequencies.start_freq_index +
- 1;
- times_and_frequencies.freq_axis = std::vector<double>(
- full_freq_axis_h5.begin() + times_and_frequencies.start_freq_index,
- full_freq_axis_h5.begin() + times_and_frequencies.start_freq_index +
- times_and_frequencies.num_freqs);
+ tf.start_freq_index = GetFreqIndex(frequencies.front());
+ tf.num_freqs = GetFreqIndex(frequencies.back()) - tf.start_freq_index + 1;
+ tf.freq_axis = std::vector<double>(
+ full_freq_axis_h5.begin() + tf.start_freq_index,
+ full_freq_axis_h5.begin() + tf.start_freq_index + tf.num_freqs);
}
if (HasAxis("pol")) {
const size_t num_pol_h5 = GetAxis("pol").size;
@@ -402,7 +408,7 @@ TimesAndFrequencies SolTab::GetTimesAndFrequencies(
}
}
- return times_and_frequencies;
+ return tf;
}
std::vector<double> SolTab::GetSubArray(
diff --git a/src/h5parm/test/tsoltab.cc b/src/h5parm/test/tsoltab.cc
index bf78cc4c39779945f5a29624edcdbe5e998f84e0..6bb1035e4fa4b9fae0f455a1198751333620a769 100644
--- a/src/h5parm/test/tsoltab.cc
+++ b/src/h5parm/test/tsoltab.cc
@@ -10,6 +10,7 @@
#include <H5Cpp.h>
using schaapcommon::h5parm::AxisInfo;
+using schaapcommon::h5parm::SelectAxisValues;
using schaapcommon::h5parm::SolTab;
namespace {
@@ -188,4 +189,141 @@ BOOST_FIXTURE_TEST_CASE(set_antennas, H5Fixture) {
BOOST_TEST(soltab.GetAntIndex(kNewNames[0]) == 0);
}
+BOOST_AUTO_TEST_CASE(select_axis_with_nearest_single_value) {
+ auto CheckSingle = [](double search_value, double expected_value,
+ size_t expected_index) {
+ size_t start_index;
+ constexpr bool kNearest = true;
+ const std::vector<double> result =
+ SelectAxisValues({4.0, 8.0, 12.0, 16.0, 20.0}, search_value,
+ search_value, kNearest, start_index);
+ BOOST_TEST(result == std::vector<double>{expected_value},
+ boost::test_tools::per_element());
+ BOOST_CHECK_EQUAL(start_index, expected_index);
+ };
+
+ // Exact matches
+ CheckSingle(4.0, 4.0, 0);
+ CheckSingle(16.0, 16.0, 3);
+ CheckSingle(20.0, 20.0, 4);
+
+ // Rounding
+ CheckSingle(5.0, 4.0, 0);
+ CheckSingle(14.1, 16.0, 3);
+ CheckSingle(17.9, 16.0, 3);
+ CheckSingle(19.0, 20.0, 4);
+
+ // Boundaries
+ CheckSingle(3.0, 4.0, 0);
+ CheckSingle(21.0, 20.0, 4);
+}
+
+BOOST_AUTO_TEST_CASE(select_axis_with_surrounding_single_value) {
+ auto CheckSingle = [](double search_value,
+ const std::vector<double>& expected_values,
+ size_t expected_index) {
+ size_t start_index;
+ constexpr bool kNearest = false;
+ const std::vector<double> result =
+ SelectAxisValues({4.0, 8.0, 12.0, 16.0, 20.0}, search_value,
+ search_value, kNearest, start_index);
+ BOOST_TEST(result == expected_values, boost::test_tools::per_element());
+ BOOST_CHECK_EQUAL(start_index, expected_index);
+ };
+
+ // Exact matches
+ CheckSingle(4.0, {4.0}, 0);
+ CheckSingle(16.0, {16.0}, 3);
+ CheckSingle(20.0, {20.0}, 4);
+
+ // Not exact matches
+ CheckSingle(5.0, {4.0, 8.0}, 0);
+ CheckSingle(14.1, {12.0, 16.0}, 2);
+ CheckSingle(17.9, {16.0, 20.0}, 3);
+ CheckSingle(19.0, {16.0, 20.0}, 3);
+
+ // Boundaries
+ CheckSingle(3.0, {4.0}, 0);
+ CheckSingle(21.0, {20.0}, 4);
+}
+
+BOOST_AUTO_TEST_CASE(select_axis_with_nearest_multiple_values) {
+ auto Check = [](double requested_start, double requested_end,
+ const std::vector<double>& expected_values,
+ size_t expected_index) {
+ size_t start_index;
+ constexpr bool kNearest = true;
+ const std::vector<double> result =
+ SelectAxisValues({4.0, 8.0, 12.0, 16.0, 20.0}, requested_start,
+ requested_end, kNearest, start_index);
+ BOOST_TEST(result == expected_values, boost::test_tools::per_element());
+ BOOST_CHECK_EQUAL(start_index, expected_index);
+ };
+
+ // Exact matches
+ Check(4.0, 8.0, {4.0, 8.0}, 0);
+ Check(8.0, 16.0, {8.0, 12.0, 16.0}, 1);
+ Check(16.0, 20.0, {16.0, 20.0}, 3);
+ Check(4.0, 20.0, {4.0, 8.0, 12.0, 16.0, 20.0}, 0);
+
+ // Rounding
+ // Both from lower
+ Check(7.0, 7.5, {8.0}, 1);
+ Check(7.0, 11.0, {8.0, 12.0}, 1);
+ // First from lower, last from higher
+ Check(7.0, 17.0, {8.0, 12.0, 16.0}, 1);
+ Check(7.0, 9.0, {8.0}, 1);
+ Check(7.0, 17.0, {8.0, 12.0, 16.0}, 1);
+ // First from higher, last from lower
+ Check(9.0, 15.0, {8.0, 12.0, 16.0}, 1);
+ Check(5.0, 17.0, {4.0, 8.0, 12.0, 16.0}, 0);
+ // First from higher, last from higher
+ Check(13.0, 17.0, {12.0, 16.0}, 2);
+ Check(5.0, 17.0, {4.0, 8.0, 12.0, 16.0}, 0);
+
+ // Boundaries
+ Check(3.0, 21.0, {4.0, 8.0, 12.0, 16.0, 20.0}, 0);
+ Check(3.0, 7.0, {4.0, 8.0}, 0);
+ Check(3.0, 4.0, {4.0}, 0);
+ Check(16.0, 21.0, {16.0, 20.0}, 3);
+ Check(17.0, 21.0, {16.0, 20.0}, 3);
+ Check(20.0, 21.0, {20.0}, 4);
+}
+
+BOOST_AUTO_TEST_CASE(select_axis_with_surrounding_multiple_values) {
+ auto Check = [](double requested_start, double requested_end,
+ const std::vector<double>& expected_values,
+ size_t expected_index) {
+ size_t start_index;
+ constexpr bool kNearest = false;
+ const std::vector<double> result =
+ SelectAxisValues({4.0, 8.0, 12.0, 16.0, 20.0}, requested_start,
+ requested_end, kNearest, start_index);
+ BOOST_TEST(result == expected_values, boost::test_tools::per_element());
+ BOOST_CHECK_EQUAL(start_index, expected_index);
+ };
+
+ // Exact matches
+ Check(4.0, 8.0, {4.0, 8.0}, 0);
+ Check(8.0, 16.0, {8.0, 12.0, 16.0}, 1);
+ Check(16.0, 20.0, {16.0, 20.0}, 3);
+ Check(4.0, 20.0, {4.0, 8.0, 12.0, 16.0, 20.0}, 0);
+
+ // Not exact matches
+ Check(7.0, 7.5, {4.0, 8.0}, 0);
+ Check(7.0, 11.0, {4.0, 8.0, 12.0}, 0);
+ Check(7.0, 17.0, {4.0, 8.0, 12.0, 16.0, 20.0}, 0);
+ Check(9.0, 15.0, {8.0, 12.0, 16.0}, 1);
+ Check(13.0, 17.0, {12.0, 16.0, 20.0}, 2);
+ Check(17.0, 19.0, {16.0, 20.0}, 3);
+
+ // Boundaries
+ Check(3.0, 21.0, {4.0, 8.0, 12.0, 16.0, 20.0}, 0);
+ Check(3.0, 7.0, {4.0, 8.0}, 0);
+ Check(3.0, 4.0, {4.0}, 0);
+ Check(16.0, 21.0, {16.0, 20.0}, 3);
+ Check(17.0, 21.0, {16.0, 20.0}, 3);
+ Check(20.0, 21.0, {20.0}, 4);
+}
+
BOOST_AUTO_TEST_SUITE_END()