Skip to content
Snippets Groups Projects
Commit 1725bcfc authored by Andre Offringa's avatar Andre Offringa
Browse files

Merge branch 'dont-store-parm-values' into 'master' 

Avoid storing duplicate buffers in JonesParameters class

See merge request !129
parents 1e402bb8 fb854e0f
No related branches found
No related tags found
1 merge request!129Avoid storing duplicate buffers in JonesParameters class
Pipeline #61371 failed
Stage: versioning
Stage: prepare
Stage: linting
Stage: build
Stage: test
Hide whitespace changes
Inline Side-by-side
include/schaapcommon/h5parm/jonesparameters.h
+ 17
− 16
View file @ 1725bcfc
......@@ -64,7 +64,11 @@ class JonesParameters {
* \param correct_type Correction type of the Jones matrix
* \param interpolation_type Interpolation type of the Jones matrix
* \param direction Direction number in the H5parm
* \param parm_values Parameter values (e.g. TEC values)
* \param parm_values Parameter values (e.g. TEC values). Inner vector
* has dimension n_times * n_frequencies, the middle vector has
* dimension n_antennas, outer vector has dimension n_parameters
* (e.g. phase and amplitude). These are the parameters as they are
* stored (e.g. TEC values).
* \param invert (optional default=false) Invert the parameters
* \param sigma_mmse (optional default=0.) Minimum mean square error parameter
* (remnant from BBS times, leave at 0 unless you know what you're doing).
......@@ -168,8 +172,10 @@ class JonesParameters {
* ONLY have an effect on RotationMeasure and RotationAngle. Other effects
* have to be inverted explicitly by calling Invert()
*/
void MakeComplex(size_t ant, const std::vector<double>& freqs,
GainType correct_type, bool invert = false);
void MakeComplex(
const std::vector<std::vector<std::vector<double>>>& parm_values,
size_t ant, const std::vector<double>& freqs, GainType correct_type,
bool invert = false);
/**
* Get the number of parameters for a given @param correct_type.
......@@ -194,13 +200,14 @@ class JonesParameters {
* \param correct_type Correction type of the Jones matrix
* \param interpolation_type Interpolation type of the Jones matrix
*/
void FillParmValues(schaapcommon::h5parm::SolTab* sol_tab,
schaapcommon::h5parm::SolTab* sol_tab2,
const std::vector<double>& freqs,
const std::vector<double>& times,
const std::vector<std::string>& antenna_names, size_t ant,
GainType gain_type, InterpolationType interpolation_type,
hsize_t direction);
void FillParmValues(
std::vector<std::vector<std::vector<double>>>& parm_values,
schaapcommon::h5parm::SolTab* sol_tab,
schaapcommon::h5parm::SolTab* sol_tab2, const std::vector<double>& freqs,
const std::vector<double>& times,
const std::vector<std::string>& antenna_names, size_t ant,
GainType gain_type, InterpolationType interpolation_type,
hsize_t direction);
/**
* Replace values by NaN on places where weight is zero
......@@ -220,12 +227,6 @@ class JonesParameters {
static void Invert(casacore::Cube<std::complex<float>>& parms,
float sigma_mmse, GainType gain_type);
/// Parameter values, inner vector has dimension num_times * num_frequencies,
/// the middle vector has dimension number_antennas, outer vector has
/// dimension num_parameters (e.g. phase and amplitude)
/// These are the paramters as they are stored (e.g. TEC values)
/// TODO: This probably does not need to be a member
std::vector<std::vector<std::vector<double>>> parm_values_;
/// Stored Jones matrices, dimensions (nparms, nantenna, ntime*nfreq),
/// frequency varies fastest. nparms is 2 for diagonal, 4 for full jones
/// parameters.
......
src/h5parm/jonesparameters.cc
+ 54
− 51
View file @ 1725bcfc
......@@ -20,12 +20,11 @@ JonesParameters::JonesParameters(
const std::vector<std::string>& antenna_names, GainType gain_type,
InterpolationType interpolation_type, hsize_t direction,
std::vector<std::vector<std::vector<double>>>&& parm_values, bool invert,
float sigma_mmse)
: parm_values_(std::move(parm_values)) {
float sigma_mmse) {
const unsigned int num_parms = GetNParms(gain_type);
parms_.resize(num_parms, antenna_names.size(), freqs.size() * times.size());
for (size_t ant = 0; ant < antenna_names.size(); ++ant) {
MakeComplex(ant, freqs, gain_type, invert);
MakeComplex(parm_values, ant, freqs, gain_type, invert);
}
if (invert) {
Invert(parms_, sigma_mmse, gain_type);
......@@ -72,15 +71,16 @@ JonesParameters::JonesParameters(
if (parm_size == 0U) {
parm_size = GetNParmValues(gain_type);
}
parm_values_.resize(parm_size);
for (auto& parm_values : parm_values_) {
std::vector<std::vector<std::vector<double>>> parm_values(parm_size);
for (auto& parm_values : parm_values) {
parm_values.resize(antenna_names.size());
}
for (size_t ant = 0; ant < antenna_names.size(); ++ant) {
try {
FillParmValues(sol_tab, sol_tab2, freqs, times, antenna_names, ant,
gain_type, interpolation_type, direction);
MakeComplex(ant, freqs, gain_type, invert);
FillParmValues(parm_values, sol_tab, sol_tab2, freqs, times,
antenna_names, ant, gain_type, interpolation_type,
direction);
MakeComplex(parm_values, ant, freqs, gain_type, invert);
} catch (const std::exception& e) {
if (std::string(e.what()).rfind("SolTab has no element", 0) != 0) {
throw;
......@@ -119,57 +119,59 @@ JonesParameters::JonesParameters(
}
}
void JonesParameters::MakeComplex(size_t ant, const std::vector<double>& freqs,
GainType gain_type, bool invert) {
void JonesParameters::MakeComplex(
const std::vector<std::vector<std::vector<double>>>& parm_values,
size_t ant, const std::vector<double>& freqs, GainType gain_type,
bool invert) {
for (unsigned int tf = 0; tf < parms_.shape()[2]; ++tf) {
const double freq = freqs[tf % freqs.size()];
switch (gain_type) {
case GainType::kDiagonalComplex:
parms_(0, ant, tf) =
casacore::polar(parm_values_[0][ant][tf], parm_values_[1][ant][tf]);
casacore::polar(parm_values[0][ant][tf], parm_values[1][ant][tf]);
parms_(1, ant, tf) =
casacore::polar(parm_values_[2][ant][tf], parm_values_[3][ant][tf]);
casacore::polar(parm_values[2][ant][tf], parm_values[3][ant][tf]);
break;
case GainType::kFullJones:
parms_(0, ant, tf) =
casacore::polar(parm_values_[0][ant][tf], parm_values_[1][ant][tf]);
casacore::polar(parm_values[0][ant][tf], parm_values[1][ant][tf]);
parms_(1, ant, tf) =
casacore::polar(parm_values_[2][ant][tf], parm_values_[3][ant][tf]);
casacore::polar(parm_values[2][ant][tf], parm_values[3][ant][tf]);
parms_(2, ant, tf) =
casacore::polar(parm_values_[4][ant][tf], parm_values_[5][ant][tf]);
casacore::polar(parm_values[4][ant][tf], parm_values[5][ant][tf]);
parms_(3, ant, tf) =
casacore::polar(parm_values_[6][ant][tf], parm_values_[7][ant][tf]);
casacore::polar(parm_values[6][ant][tf], parm_values[7][ant][tf]);
break;
case GainType::kScalarComplex:
parms_(0, ant, tf) =
casacore::polar(parm_values_[0][ant][tf], parm_values_[1][ant][tf]);
casacore::polar(parm_values[0][ant][tf], parm_values[1][ant][tf]);
parms_(1, ant, tf) = parms_(0, ant, tf);
break;
case GainType::kTec:
parms_(0, ant, tf) = casacore::polar(
1., parm_values_[0][ant][tf] * -8.44797245e9 / freq);
if (parm_values_.size() == 1) { // No TEC:0, only TEC:
parms_(0, ant, tf) =
casacore::polar(1., parm_values[0][ant][tf] * -8.44797245e9 / freq);
if (parm_values.size() == 1) { // No TEC:0, only TEC:
parms_(1, ant, tf) = casacore::polar(
1., parm_values_[0][ant][tf] * -8.44797245e9 / freq);
1., parm_values[0][ant][tf] * -8.44797245e9 / freq);
} else { // TEC:0 and TEC:1
parms_(1, ant, tf) = casacore::polar(
1., parm_values_[1][ant][tf] * -8.44797245e9 / freq);
1., parm_values[1][ant][tf] * -8.44797245e9 / freq);
}
break;
case GainType::kClock:
parms_(0, ant, tf) = casacore::polar(
1., parm_values_[0][ant][tf] * freq * casacore::C::_2pi);
if (parm_values_.size() == 1) { // No Clock:0, only Clock:
1., parm_values[0][ant][tf] * freq * casacore::C::_2pi);
if (parm_values.size() == 1) { // No Clock:0, only Clock:
parms_(1, ant, tf) = casacore::polar(
1., parm_values_[0][ant][tf] * freq * casacore::C::_2pi);
1., parm_values[0][ant][tf] * freq * casacore::C::_2pi);
} else { // Clock:0 and Clock:1
parms_(1, ant, tf) = casacore::polar(
1., parm_values_[1][ant][tf] * freq * casacore::C::_2pi);
1., parm_values[1][ant][tf] * freq * casacore::C::_2pi);
}
break;
case GainType::kRotationAngle: {
double phi = parm_values_[0][ant][tf];
double phi = parm_values[0][ant][tf];
if (invert) {
phi = -phi;
}
......@@ -183,7 +185,7 @@ void JonesParameters::MakeComplex(size_t ant, const std::vector<double>& freqs,
case GainType::kRotationMeasure: {
const double lambda2 =
(casacore::C::c / freq) * (casacore::C::c / freq);
double chi = parm_values_[0][ant][tf] * lambda2;
double chi = parm_values[0][ant][tf] * lambda2;
if (invert) {
chi = -chi;
}
......@@ -196,38 +198,38 @@ void JonesParameters::MakeComplex(size_t ant, const std::vector<double>& freqs,
} break;
case GainType::kDiagonalPhase:
case GainType::kScalarPhase:
parms_(0, ant, tf) = casacore::polar(1., parm_values_[0][ant][tf]);
parms_(0, ant, tf) = casacore::polar(1., parm_values[0][ant][tf]);
if (gain_type == GainType::kScalarPhase) { // Same value for x and y
parms_(1, ant, tf) = casacore::polar(1., parm_values_[0][ant][tf]);
parms_(1, ant, tf) = casacore::polar(1., parm_values[0][ant][tf]);
} else { // Different value for x and y
parms_(1, ant, tf) = casacore::polar(1., parm_values_[1][ant][tf]);
parms_(1, ant, tf) = casacore::polar(1., parm_values[1][ant][tf]);
}
break;
case GainType::kDiagonalAmplitude:
case GainType::kScalarAmplitude:
parms_(0, ant, tf) = parm_values_[0][ant][tf];
parms_(0, ant, tf) = parm_values[0][ant][tf];
if (gain_type ==
GainType::kScalarAmplitude) { // Same value for x and y
parms_(1, ant, tf) = parm_values_[0][ant][tf];
parms_(1, ant, tf) = parm_values[0][ant][tf];
} else { // Different value for x and y
parms_(1, ant, tf) = parm_values_[1][ant][tf];
parms_(1, ant, tf) = parm_values[1][ant][tf];
}
break;
case GainType::kDiagonalRealImaginary:
parms_(0, ant, tf) = std::complex<float>(parm_values_[0][ant][tf],
parm_values_[1][ant][tf]);
parms_(1, ant, tf) = std::complex<float>(parm_values_[2][ant][tf],
parm_values_[3][ant][tf]);
parms_(0, ant, tf) = std::complex<float>(parm_values[0][ant][tf],
parm_values[1][ant][tf]);
parms_(1, ant, tf) = std::complex<float>(parm_values[2][ant][tf],
parm_values[3][ant][tf]);
break;
case GainType::kFullJonesRealImaginary:
parms_(0, ant, tf) = std::complex<float>(parm_values_[0][ant][tf],
parm_values_[1][ant][tf]);
parms_(1, ant, tf) = std::complex<float>(parm_values_[2][ant][tf],
parm_values_[3][ant][tf]);
parms_(2, ant, tf) = std::complex<float>(parm_values_[4][ant][tf],
parm_values_[5][ant][tf]);
parms_(3, ant, tf) = std::complex<float>(parm_values_[6][ant][tf],
parm_values_[7][ant][tf]);
parms_(0, ant, tf) = std::complex<float>(parm_values[0][ant][tf],
parm_values[1][ant][tf]);
parms_(1, ant, tf) = std::complex<float>(parm_values[2][ant][tf],
parm_values[3][ant][tf]);
parms_(2, ant, tf) = std::complex<float>(parm_values[4][ant][tf],
parm_values[5][ant][tf]);
parms_(3, ant, tf) = std::complex<float>(parm_values[6][ant][tf],
parm_values[7][ant][tf]);
}
}
}
......@@ -270,6 +272,7 @@ unsigned int JonesParameters::GetNParmValues(GainType gain_type) {
}
void JonesParameters::FillParmValues(
std::vector<std::vector<std::vector<double>>>& parm_values,
schaapcommon::h5parm::SolTab* sol_tab,
schaapcommon::h5parm::SolTab* sol_tab2, const std::vector<double>& freqs,
const std::vector<double>& times,
......@@ -294,20 +297,20 @@ void JonesParameters::FillParmValues(
if (gain_type == GainType::kFullJones ||
gain_type == GainType::kDiagonalComplex ||
gain_type == GainType::kScalarComplex) {
for (size_t pol = 0; pol < parm_values_.size() / 2; ++pol) {
for (size_t pol = 0; pol < parm_values.size() / 2; ++pol) {
// Place amplitude in even and phase in odd elements
parm_values_[pol * 2][ant] = get_flagged_values(sol_tab, ant_name, pol);
parm_values[pol * 2][ant] = get_flagged_values(sol_tab, ant_name, pol);
if (sol_tab2 == nullptr) {
throw std::runtime_error(
"soltab2 cannot be a nullpointer for correct_type=FULLJONES and "
"correct_type=GAIN");
}
parm_values_[pol * 2 + 1][ant] =
parm_values[pol * 2 + 1][ant] =
get_flagged_values(sol_tab2, ant_name, pol);
}
} else {
for (size_t pol = 0; pol < parm_values_.size(); ++pol) {
parm_values_[pol][ant] = get_flagged_values(sol_tab, ant_name, pol);
for (size_t pol = 0; pol < parm_values.size(); ++pol) {
parm_values[pol][ant] = get_flagged_values(sol_tab, ant_name, pol);
}
}
}
......
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