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

Avoid storing duplicate buffers in JonesParameters class

parent 1e402bb8
No related branches found
No related tags found
1 merge request!129Avoid storing duplicate buffers in JonesParameters class
Hide whitespace changes
Inline Side-by-side
include/schaapcommon/h5parm/jonesparameters.h
+ 17
− 16
View file @ fb854e0f
......@@ -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 @ fb854e0f
......@@ -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