Skip to content
Snippets Groups Projects
Commit e00188db authored by Tammo Jan Dijkema's avatar Tammo Jan Dijkema
Browse files

Merge branch 'more_extrapolate2' into 'master' 

Add some more extrapolation in h5parm freqs

See merge request !48
parents 99df042f 87f21e94
No related branches found
No related tags found
1 merge request!48Add some more extrapolation in h5parm freqs
Pipeline #13493 passed
Stage: prepare
Stage: build
Stage: test
Hide whitespace changes
Inline Side-by-side
include/schaapcommon/h5parm/soltab.h
+ 4
− 2
View file @ e00188db
......@@ -75,8 +75,10 @@ class SolTab : private H5::Group {
/// Get the values of a string-valued axis (e.g. "dir" or "pol")
std::vector<std::string> GetStringAxis(const std::string& axis_name);
/// Get the index of freq, using nearest neighbor
/// This assumes that the frequencies are in increasing order.
/// Get the index of freq, using nearest neighbor.
/// This function assumes that the frequencies are in increasing order.
/// @throw std::runtime_error If the frequency is less than a full frequency
/// width below the lowest or above the highest frequency
hsize_t GetFreqIndex(double freq) const;
/// Get the index of a time. Matches with 0.5*timeInterval
......
src/h5parm/soltab.cc
+ 12
− 10
View file @ e00188db
......@@ -99,7 +99,7 @@ void SolTab::SetValues(const std::vector<double>& vals,
const std::string& history) {
// Convert axes to comma separated string, fill dims
size_t expectedsize = 1;
std::string axesstr = axes_[0].name;
std::string axesstr = axes_.front().name;
std::vector<hsize_t> dims(axes_.size());
for (unsigned int i = 0; i < axes_.size(); ++i) {
dims[i] = axes_[i].size;
......@@ -276,8 +276,8 @@ std::vector<double> SolTab::GetValuesOrWeights(const std::string& val_or_weight,
}
if (HasAxis("freq")) {
std::vector<double> full_freq_axis_h5 = GetRealAxis("freq");
freq_start = GetFreqIndex(freqs[0]);
num_freq_h5 = GetFreqIndex(freqs[freqs.size() - 1]) - freq_start + 1;
freq_start = GetFreqIndex(freqs.front());
num_freq_h5 = GetFreqIndex(freqs.back()) - freq_start + 1;
freq_axis_h5 = std::vector<double>(
full_freq_axis_h5.begin() + freq_start,
full_freq_axis_h5.begin() + freq_start + num_freq_h5);
......@@ -291,7 +291,6 @@ std::vector<double> SolTab::GetValuesOrWeights(const std::string& val_or_weight,
" polarizations are in there.");
}
}
std::vector<double> h5values =
GetValuesOrWeights(val_or_weight, ant_name, start_time_slot, num_time_h5,
1, freq_start, num_freq_h5, 1, pol, dir);
......@@ -490,16 +489,19 @@ hsize_t SolTab::GetFreqIndex(double freq) const {
return 0;
}
std::vector<double> freqs = GetRealAxis("freq");
double freq_interval = GetFreqInterval(0);
// Half a cell width before the first frequency
if (abs(freqs[0] - freq) < 0.501 * freq_interval) {
// A full cell width before the first frequency
if (freq < freqs.front() - freq_interval) {
throw std::runtime_error("Frequency " + std::to_string(freq) +
" not found in " + GetName());
}
if (freq < freqs.front()) {
return 0;
}
// No assumptions on regular spacing here
for (size_t i = 0; i < freqs.size() - 1; ++i) {
if (freqs[i] - 0.001 <= freq && freq < freqs[i + 1]) { // Some tolerance
if (freq < freqs[i + 1]) {
// Nearest neighbor: i or i+1
if (freq - freqs[i] < freqs[i + 1] - freq) {
return i;
......@@ -509,9 +511,9 @@ hsize_t SolTab::GetFreqIndex(double freq) const {
}
}
// Half a cell width after the last frequency
// A full cell width after the last frequency
freq_interval = GetFreqInterval(freqs.size() - 2);
if (abs(freqs[freqs.size() - 1] - freq) < 0.501 * freq_interval) {
if (freq < freqs.back() + freq_interval) {
return freqs.size() - 1;
}
......
src/h5parm/test/th5parm.cc
+ 50
− 10
View file @ e00188db
......@@ -20,6 +20,11 @@ using std::vector;
BOOST_AUTO_TEST_SUITE(h5parm)
namespace {
const size_t kNumAntennas = 3;
const size_t kNumFrequencies = 4;
void CheckAxes(SolTab& soltab, size_t ntimes) {
BOOST_CHECK_EQUAL(soltab.NumAxes(), size_t{3});
BOOST_CHECK(soltab.HasAxis("ant"));
......@@ -50,16 +55,21 @@ void InitializeH5(H5Parm& h5parm, size_t ntimes) {
axes.push_back(AxisInfo("time", ntimes));
axes.push_back(AxisInfo("bla", 1));
SolTab a = h5parm.CreateSolTab("mysol", "mytype", axes);
h5parm.CreateSolTab("mysol", "mytype", axes);
vector<AxisInfo> axes_freq;
axes_freq.push_back(AxisInfo("ant", kNumAntennas));
axes_freq.push_back(AxisInfo("freq", kNumFrequencies));
h5parm.CreateSolTab("mysolwithfreq", "mytype", axes_freq);
}
void FillData(H5Parm& h5parm, size_t ntimes) {
SolTab soltab = h5parm.GetSolTab("mysol");
// Add some data
vector<double> vals(3 * ntimes);
vector<double> weights(3 * ntimes);
for (size_t ant = 0; ant < 3; ++ant) {
vector<double> vals(kNumAntennas * ntimes);
vector<double> weights(kNumAntennas * ntimes);
for (size_t ant = 0; ant < kNumAntennas; ++ant) {
for (size_t time = 0; time < ntimes; ++time) {
vals[ant * ntimes + time] = 10 * ant + time;
weights[ant * ntimes + time] = 0.4;
......@@ -81,11 +91,27 @@ void FillData(H5Parm& h5parm, size_t ntimes) {
times.push_back(57878.5 + 2.0 * time);
}
soltab.SetTimes(times);
}
void FillDataWithFreqAxis(H5Parm& h5parm) {
SolTab soltab = h5parm.GetSolTab("mysolwithfreq");
// Add some data
vector<double> vals(kNumAntennas * kNumFrequencies, 1.0);
vector<double> weights(kNumAntennas * kNumFrequencies, 1.0);
soltab.SetValues(vals, weights, "CREATE with SchaapCommon-test");
// Add metadata for freqs;
vector<double> freqs{130e6, 131e6, 135e6, 137e6};
soltab.SetFreqs(freqs);
// Add metadata for stations
const std::vector<std::string> someAntNames{"Antenna1", "Antenna12",
"Antenna42"};
soltab.SetAntennas(someAntNames);
}
} // namespace
BOOST_AUTO_TEST_CASE(create) {
constexpr size_t ntimes = 7;
......@@ -102,7 +128,7 @@ BOOST_AUTO_TEST_CASE(create) {
InitializeH5(h5parm, ntimes);
// Check that the soltab exists
BOOST_CHECK_EQUAL(h5parm.NumSolTabs(), size_t{1});
BOOST_CHECK_EQUAL(h5parm.NumSolTabs(), size_t{2});
BOOST_CHECK(h5parm.HasSolTab("mysol"));
// Check the axes
......@@ -116,6 +142,7 @@ struct H5Fixture {
H5Parm h5parm("tH5Parm_tmp.h5", true);
InitializeH5(h5parm, 7);
FillData(h5parm, 7);
FillDataWithFreqAxis(h5parm);
}
~H5Fixture() { remove("tH5Parm_tmp.h5"); }
......@@ -130,7 +157,7 @@ BOOST_FIXTURE_TEST_CASE(open, H5Fixture) {
constexpr size_t ntimes = 7;
H5Parm h5parm("tH5Parm_tmp.h5", false, false, "sol000");
BOOST_CHECK_EQUAL(h5parm.GetSolSetName(), "sol000");
BOOST_CHECK_EQUAL(h5parm.NumSolTabs(), size_t{1});
BOOST_CHECK_EQUAL(h5parm.NumSolTabs(), size_t{2});
BOOST_CHECK(h5parm.HasSolTab("mysol"));
BOOST_CHECK(!h5parm.HasSolTab("nonexistingsol"));
......@@ -164,10 +191,7 @@ BOOST_FIXTURE_TEST_CASE(open, H5Fixture) {
BOOST_CHECK_EQUAL(antennas[1], "Antenna12");
BOOST_CHECK_EQUAL(antennas[2], "Antenna123");
BOOST_CHECK_CLOSE(soltab.GetFreqInterval(0), 1e6, 1e-8);
BOOST_CHECK_CLOSE(soltab.GetFreqInterval(1), 4e6, 1e-8);
BOOST_CHECK_CLOSE(soltab.GetFreqInterval(2), 2e6, 1e-8);
// Check grid interpolation
vector<double> freqs{130e6, 131e6};
vector<double> times;
......@@ -200,6 +224,22 @@ BOOST_FIXTURE_TEST_CASE(open, H5Fixture) {
1e-8);
}
}
soltab = h5parm.GetSolTab("mysolwithfreq");
BOOST_CHECK_CLOSE(soltab.GetFreqInterval(0), 1.0e6, 1.0e-8);
BOOST_CHECK_CLOSE(soltab.GetFreqInterval(1), 4.0e6, 1.0e-8);
BOOST_CHECK_CLOSE(soltab.GetFreqInterval(2), 2.0e6, 1.0e-8);
BOOST_CHECK_THROW(soltab.GetFreqIndex(128.0e6),
std::runtime_error); // Too far from lowest
BOOST_CHECK_EQUAL(soltab.GetFreqIndex(129.1e6), 0); // closest to 130e6
BOOST_CHECK_EQUAL(soltab.GetFreqIndex(130.4e6), 0); // closest to 130e6
BOOST_CHECK_EQUAL(soltab.GetFreqIndex(130.6e6), 1); // closest to 131e6
BOOST_CHECK_EQUAL(soltab.GetFreqIndex(136.1e6), 3); // closest to 137e6
BOOST_CHECK_EQUAL(soltab.GetFreqIndex(137.0e6), 3); // closest to 137e6
BOOST_CHECK_EQUAL(soltab.GetFreqIndex(137.8e6), 3); // closest to 137e6
BOOST_CHECK_THROW(soltab.GetFreqIndex(150.0e6),
std::runtime_error); // Too far from highest
}
BOOST_AUTO_TEST_SUITE_END()
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