diff --git a/base/Apply.cc b/base/Apply.cc
index 4b08fa9f0dcfd94082c1de82ab91293ce7018c80..4ad026b79e932c6297cb6733e54ac5c50925ebbf 100644
--- a/base/Apply.cc
+++ b/base/Apply.cc
@@ -11,7 +11,7 @@ namespace dp3 {
namespace base {
void apply(size_t nBaseline, size_t nChannel, const_cursor<Baseline> baselines,
- const_cursor<double> coeff, cursor<dcomplex> data) {
+ const_cursor<double> coeff, cursor<std::complex<double>> data) {
for (size_t bl = 0; bl < nBaseline; ++bl) {
const size_t p = baselines->first;
const size_t q = baselines->second;
@@ -19,36 +19,36 @@ void apply(size_t nBaseline, size_t nChannel, const_cursor<Baseline> baselines,
if (p != q) {
// Jones matrix for station P.
coeff.forward(1, p);
- const dcomplex Jp_00(coeff[0], coeff[1]);
- const dcomplex Jp_01(coeff[2], coeff[3]);
- const dcomplex Jp_10(coeff[4], coeff[5]);
- const dcomplex Jp_11(coeff[6], coeff[7]);
+ const std::complex<double> Jp_00(coeff[0], coeff[1]);
+ const std::complex<double> Jp_01(coeff[2], coeff[3]);
+ const std::complex<double> Jp_10(coeff[4], coeff[5]);
+ const std::complex<double> Jp_11(coeff[6], coeff[7]);
coeff.backward(1, p);
// Jones matrix for station Q, conjugated.
coeff.forward(1, q);
- const dcomplex Jq_00(coeff[0], -coeff[1]);
- const dcomplex Jq_01(coeff[2], -coeff[3]);
- const dcomplex Jq_10(coeff[4], -coeff[5]);
- const dcomplex Jq_11(coeff[6], -coeff[7]);
+ const std::complex<double> Jq_00(coeff[0], -coeff[1]);
+ const std::complex<double> Jq_01(coeff[2], -coeff[3]);
+ const std::complex<double> Jq_10(coeff[4], -coeff[5]);
+ const std::complex<double> Jq_11(coeff[6], -coeff[7]);
coeff.backward(1, q);
// Compute (Jp x conj(Jq)) * vec(data), where 'x' denotes the
// Kronecker product.
for (size_t ch = 0; ch < nChannel; ++ch) {
// Fetch visibilities.
- const dcomplex xx = data[0];
- const dcomplex xy = data[1];
- const dcomplex yx = data[2];
- const dcomplex yy = data[3];
+ const std::complex<double> xx = data[0];
+ const std::complex<double> xy = data[1];
+ const std::complex<double> yx = data[2];
+ const std::complex<double> yy = data[3];
// Precompute terms involving conj(Jq) and data. Each term
// appears twice in the computation of (Jp x conj(Jq))
// * vec(data).
- const dcomplex Jq_00xx_01xy = Jq_00 * xx + Jq_01 * xy;
- const dcomplex Jq_00yx_01yy = Jq_00 * yx + Jq_01 * yy;
- const dcomplex Jq_10xx_11xy = Jq_10 * xx + Jq_11 * xy;
- const dcomplex Jq_10yx_11yy = Jq_10 * yx + Jq_11 * yy;
+ const std::complex<double> Jq_00xx_01xy = Jq_00 * xx + Jq_01 * xy;
+ const std::complex<double> Jq_00yx_01yy = Jq_00 * yx + Jq_01 * yy;
+ const std::complex<double> Jq_10xx_11xy = Jq_10 * xx + Jq_11 * xy;
+ const std::complex<double> Jq_10yx_11yy = Jq_10 * yx + Jq_11 * yy;
// Compute (Jp x conj(Jq)) * vec(data) from the precomputed
// terms.
diff --git a/base/Cursor.h b/base/Cursor.h
index 2571ed1de743e021c953257ae5515720e7de376d..4c1221ac5e6c6167ce7ccb631ae6b4dffe24639a 100644
--- a/base/Cursor.h
+++ b/base/Cursor.h
@@ -3,15 +3,12 @@
// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
-#ifndef DPPP_CURSOR_H
-#define DPPP_CURSOR_H
+#ifndef DP3_BASE_CURSOR_H_
+#define DP3_BASE_CURSOR_H_
#include <cassert>
#include <complex>
-typedef std::complex<double> dcomplex;
-typedef std::complex<float> fcomplex;
-
namespace dp3 {
namespace base {
diff --git a/base/EstimateMixed.cc b/base/EstimateMixed.cc
index 99ff743cf97d0d90d3070e1dda7b67798959c685..597f733d38abc1da98b12cfa7237ef8d29304c20 100644
--- a/base/EstimateMixed.cc
+++ b/base/EstimateMixed.cc
@@ -49,10 +49,11 @@ void makeIndex(std::size_t n_direction, std::size_t n_station,
bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
size_t nChannel, const_cursor<Baseline> baselines,
- vector<const_cursor<fcomplex>> data,
- vector<const_cursor<dcomplex>> model, const_cursor<bool> flag,
- const_cursor<float> weight, const_cursor<dcomplex> mix,
- double *unknowns, size_t maxiter) {
+ std::vector<const_cursor<std::complex<float>>> data,
+ std::vector<const_cursor<std::complex<double>>> model,
+ const_cursor<bool> flag, const_cursor<float> weight,
+ const_cursor<std::complex<double>> mix, double *unknowns,
+ size_t maxiter) {
assert(data.size() == nDirection && model.size() == nDirection);
// Initialize LSQ solver.
@@ -70,8 +71,8 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
vector<unsigned int> dIndex(4 * nPartial);
// Allocate space for intermediate results.
- vector<dcomplex> M(nDirection * 4), dM(nDirection * 16);
- vector<double> dR(nPartial), dI(nPartial);
+ std::vector<std::complex<double>> M(nDirection * 4), dM(nDirection * 16);
+ std::vector<double> dR(nPartial), dI(nPartial);
// Iterate until convergence.
size_t nIterations = 0;
@@ -88,30 +89,30 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
for (size_t dr = 0; dr < nDirection; ++dr) {
// Jones matrix for station P.
const double *Jp = &(unknowns[dr * nStation * 8 + p * 8]);
- const dcomplex Jp_00(Jp[0], Jp[1]);
- const dcomplex Jp_01(Jp[2], Jp[3]);
- const dcomplex Jp_10(Jp[4], Jp[5]);
- const dcomplex Jp_11(Jp[6], Jp[7]);
+ const std::complex<double> Jp_00(Jp[0], Jp[1]);
+ const std::complex<double> Jp_01(Jp[2], Jp[3]);
+ const std::complex<double> Jp_10(Jp[4], Jp[5]);
+ const std::complex<double> Jp_11(Jp[6], Jp[7]);
// Jones matrix for station Q, conjugated.
const double *Jq = &(unknowns[dr * nStation * 8 + q * 8]);
- const dcomplex Jq_00(Jq[0], -Jq[1]);
- const dcomplex Jq_01(Jq[2], -Jq[3]);
- const dcomplex Jq_10(Jq[4], -Jq[5]);
- const dcomplex Jq_11(Jq[6], -Jq[7]);
+ const std::complex<double> Jq_00(Jq[0], -Jq[1]);
+ const std::complex<double> Jq_01(Jq[2], -Jq[3]);
+ const std::complex<double> Jq_10(Jq[4], -Jq[5]);
+ const std::complex<double> Jq_11(Jq[6], -Jq[7]);
// Fetch model visibilities for the current direction.
- const dcomplex xx = model[dr][0];
- const dcomplex xy = model[dr][1];
- const dcomplex yx = model[dr][2];
- const dcomplex yy = model[dr][3];
+ const std::complex<double> xx = model[dr][0];
+ const std::complex<double> xy = model[dr][1];
+ const std::complex<double> yx = model[dr][2];
+ const std::complex<double> yy = model[dr][3];
// Precompute terms involving conj(Jq) and the model
// visibilities.
- const dcomplex Jq_00xx_01xy = Jq_00 * xx + Jq_01 * xy;
- const dcomplex Jq_00yx_01yy = Jq_00 * yx + Jq_01 * yy;
- const dcomplex Jq_10xx_11xy = Jq_10 * xx + Jq_11 * xy;
- const dcomplex Jq_10yx_11yy = Jq_10 * yx + Jq_11 * yy;
+ const std::complex<double> Jq_00xx_01xy = Jq_00 * xx + Jq_01 * xy;
+ const std::complex<double> Jq_00yx_01yy = Jq_00 * yx + Jq_01 * yy;
+ const std::complex<double> Jq_10xx_11xy = Jq_10 * xx + Jq_11 * xy;
+ const std::complex<double> Jq_10yx_11yy = Jq_10 * yx + Jq_11 * yy;
// Precompute (Jp x conj(Jq)) * vec(data), where 'x'
// denotes the Kronecker product. This is the model
@@ -152,16 +153,16 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
{
if (!flag[cr]) {
for (size_t tg = 0; tg < nDirection; ++tg) {
- dcomplex visibility(0.0, 0.0);
+ std::complex<double> visibility(0.0, 0.0);
for (size_t dr = 0; dr < nDirection; ++dr) {
// Look-up mixing weight.
- const dcomplex mix_weight = *mix;
+ const std::complex<double> mix_weight = *mix;
// Weight model visibility.
visibility += mix_weight * M[dr * 4 + cr];
// Compute weighted partial derivatives.
- dcomplex derivative(0.0, 0.0);
+ std::complex<double> derivative(0.0, 0.0);
derivative = mix_weight * dM[dr * 16 + cr * 4];
dR[dr * 8] = real(derivative); // for cr==0: Re(d/dRe(p_00)))
dI[dr * 8] = imag(derivative); // for cr==0: Re(d/dIm(p_00)))
@@ -205,16 +206,14 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
} // Source directions.
// Compute the residual.
- dcomplex residual =
- static_cast<dcomplex>(data[tg][cr]) - visibility;
+ std::complex<double> residual =
+ std::complex<double>{data[tg][cr]} - visibility;
// Update the normal equations.
solver.makeNorm(nPartial, &(dIndex[cr * nPartial]), &(dR[0]),
- static_cast<double>(weight[cr]),
- real(residual));
+ double{weight[cr]}, residual.real());
solver.makeNorm(nPartial, &(dIndex[cr * nPartial]), &(dI[0]),
- static_cast<double>(weight[cr]),
- imag(residual));
+ double{weight[cr]}, residual.imag());
// Move to next target direction.
mix.backward(1, nDirection);
diff --git a/base/EstimateMixed.h b/base/EstimateMixed.h
index fac09ab6233b77e29185d78526e38289860d165e..9dc1fb0a2e6ff1c093499c595273d7391d695862 100644
--- a/base/EstimateMixed.h
+++ b/base/EstimateMixed.h
@@ -57,10 +57,10 @@ namespace base {
/// A pointer to a buffer of unknowns of size nDirection * nStation * 8.
bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
size_t nChannel, const_cursor<Baseline> baselines,
- std::vector<const_cursor<fcomplex>> data,
- std::vector<const_cursor<dcomplex>> model,
+ std::vector<const_cursor<std::complex<float>>> data,
+ std::vector<const_cursor<std::complex<double>>> model,
const_cursor<bool> flag, const_cursor<float> weight,
- const_cursor<dcomplex> mix, double* unknowns,
+ const_cursor<std::complex<double>> mix, double* unknowns,
size_t maxiter = 50);
#ifdef HAVE_LIBDIRAC
@@ -73,10 +73,10 @@ bool estimate(size_t nDirection, size_t nStation, size_t nBaseline,
bool estimate(std::size_t n_direction, std::size_t n_station,
std::size_t n_baseline, std::size_t n_channel,
const_cursor<Baseline> baselines,
- std::vector<const_cursor<fcomplex>> data,
- std::vector<const_cursor<dcomplex>> model,
+ std::vector<const_cursor<std::complex<float>>> data,
+ std::vector<const_cursor<std::complex<double>>> model,
const_cursor<bool> flag, const_cursor<float> weight,
- const_cursor<dcomplex> mix, double* unknowns,
+ const_cursor<std::complex<double>> mix, double* unknowns,
std::size_t lbfgs_mem, double robust_nu,
std::size_t max_iter = 50);
#endif /* HAVE_LIBDIRAC */
diff --git a/base/EstimateMixedLBFGS.cc b/base/EstimateMixedLBFGS.cc
index d6112f3400c58fa621d6cb8e5a84fc21032a0e09..25efdb5dfebce2b4795f8e4660f4aede3cfe213b 100644
--- a/base/EstimateMixedLBFGS.cc
+++ b/base/EstimateMixedLBFGS.cc
@@ -18,6 +18,9 @@
#include <iostream>
#ifdef HAVE_LIBDIRAC
#include <Dirac.h>
+// Dirac.h incorrectly defines "complex".
+// -> Undefining it avoids conflicts with std::complex.
+#undef complex
#endif /* HAVE_LIBDIRAC */
using dp3::common::operator<<;
@@ -33,18 +36,18 @@ struct LBFGSData {
std::size_t n_baseline;
std::size_t n_channel;
const_cursor<Baseline> baselines;
- std::vector<const_cursor<fcomplex>> data;
- std::vector<const_cursor<dcomplex>> model;
+ std::vector<const_cursor<std::complex<float>>> data;
+ std::vector<const_cursor<std::complex<double>>> model;
const_cursor<bool> flag;
const_cursor<float> weight;
- const_cursor<dcomplex> mix;
+ const_cursor<std::complex<double>> mix;
double robust_nu;
LBFGSData(std::size_t n_dir_, std::size_t n_st_, std::size_t n_base_,
std::size_t n_chan_, const_cursor<Baseline> baselines_,
- std::vector<const_cursor<fcomplex>> data_,
- std::vector<const_cursor<dcomplex>> model_,
+ std::vector<const_cursor<std::complex<float>>> data_,
+ std::vector<const_cursor<std::complex<double>>> model_,
const_cursor<bool> flag_, const_cursor<float> weight_,
- const_cursor<dcomplex> mix_, double robust_nu_ = 2.0)
+ const_cursor<std::complex<double>> mix_, double robust_nu_ = 2.0)
: n_direction(n_dir_),
n_station(n_st_),
n_baseline(n_base_),
@@ -69,7 +72,7 @@ double cost_func(double *unknowns, int m, void *adata) {
assert(static_cast<std::size_t>(m) == t->n_direction * t->n_station * 4 * 2);
// Allocate space for intermediate results.
- std::vector<dcomplex> M(t->n_direction * 4);
+ std::vector<std::complex<double>> M(t->n_direction * 4);
double fcost = 0.0;
@@ -83,30 +86,30 @@ double cost_func(double *unknowns, int m, void *adata) {
for (std::size_t dr = 0; dr < t->n_direction; ++dr) {
// Jones matrix for station P.
const double *Jp = &(unknowns[dr * t->n_station * 8 + p * 8]);
- const dcomplex Jp_00(Jp[0], Jp[1]);
- const dcomplex Jp_01(Jp[2], Jp[3]);
- const dcomplex Jp_10(Jp[4], Jp[5]);
- const dcomplex Jp_11(Jp[6], Jp[7]);
+ const std::complex<double> Jp_00(Jp[0], Jp[1]);
+ const std::complex<double> Jp_01(Jp[2], Jp[3]);
+ const std::complex<double> Jp_10(Jp[4], Jp[5]);
+ const std::complex<double> Jp_11(Jp[6], Jp[7]);
// Jones matrix for station Q, conjugated.
const double *Jq = &(unknowns[dr * t->n_station * 8 + q * 8]);
- const dcomplex Jq_00(Jq[0], -Jq[1]);
- const dcomplex Jq_01(Jq[2], -Jq[3]);
- const dcomplex Jq_10(Jq[4], -Jq[5]);
- const dcomplex Jq_11(Jq[6], -Jq[7]);
+ const std::complex<double> Jq_00(Jq[0], -Jq[1]);
+ const std::complex<double> Jq_01(Jq[2], -Jq[3]);
+ const std::complex<double> Jq_10(Jq[4], -Jq[5]);
+ const std::complex<double> Jq_11(Jq[6], -Jq[7]);
// Fetch model visibilities for the current direction.
- const dcomplex xx = t->model[dr][0];
- const dcomplex xy = t->model[dr][1];
- const dcomplex yx = t->model[dr][2];
- const dcomplex yy = t->model[dr][3];
+ const std::complex<double> xx = t->model[dr][0];
+ const std::complex<double> xy = t->model[dr][1];
+ const std::complex<double> yx = t->model[dr][2];
+ const std::complex<double> yy = t->model[dr][3];
// Precompute terms involving conj(Jq) and the model
// visibilities.
- const dcomplex Jq_00xx_01xy = Jq_00 * xx + Jq_01 * xy;
- const dcomplex Jq_00yx_01yy = Jq_00 * yx + Jq_01 * yy;
- const dcomplex Jq_10xx_11xy = Jq_10 * xx + Jq_11 * xy;
- const dcomplex Jq_10yx_11yy = Jq_10 * yx + Jq_11 * yy;
+ const std::complex<double> Jq_00xx_01xy = Jq_00 * xx + Jq_01 * xy;
+ const std::complex<double> Jq_00yx_01yy = Jq_00 * yx + Jq_01 * yy;
+ const std::complex<double> Jq_10xx_11xy = Jq_10 * xx + Jq_11 * xy;
+ const std::complex<double> Jq_10yx_11yy = Jq_10 * yx + Jq_11 * yy;
// Precompute (Jp x conj(Jq)) * vec(data), where 'x'
// denotes the Kronecker product. This is the model
@@ -124,10 +127,10 @@ double cost_func(double *unknowns, int m, void *adata) {
if (!(t->flag[cr])) {
const double mwt = static_cast<double>(t->weight[cr]);
for (std::size_t tg = 0; tg < t->n_direction; ++tg) {
- dcomplex visibility(0.0, 0.0);
+ std::complex<double> visibility(0.0, 0.0);
for (std::size_t dr = 0; dr < t->n_direction; ++dr) {
// Look-up mixing weight.
- const dcomplex mix_weight = *(t->mix);
+ const std::complex<double> mix_weight = *(t->mix);
// Weight model visibility.
visibility += mix_weight * M[dr * 4 + cr];
@@ -137,17 +140,17 @@ double cost_func(double *unknowns, int m, void *adata) {
} // Source directions.
// Compute the residual.
- dcomplex residual =
- static_cast<dcomplex>(t->data[tg][cr]) - visibility;
+ std::complex<double> residual =
+ std::complex<double>{t->data[tg][cr]} - visibility;
// sum up cost
// For reference: Gaussian cost is
- // fcost += mwt * (real(residual) * real(residual) +
- // imag(residual) * imag(residual));
+ // fcost += mwt * (residual.real() * residual.real() +
+ // residual.imag() * residual.imag());
// Robust cost function
- fcost += std::log(1.0 + mwt * real(residual) * real(residual) /
+ fcost += std::log(1.0 + mwt * residual.real() * residual.real() /
t->robust_nu);
- fcost += std::log(1.0 + mwt * imag(residual) * imag(residual) /
+ fcost += std::log(1.0 + mwt * residual.imag() * residual.imag() /
t->robust_nu);
// Move to next target direction.
@@ -236,8 +239,8 @@ void grad_func(double *unknowns, double *grad, int m, void *adata) {
t->n_direction * 8; // note: this for real,imag separately
// Allocate space for intermediate results.
- std::vector<dcomplex> M(t->n_direction * 4);
- std::vector<dcomplex> dM(t->n_direction * 16);
+ std::vector<std::complex<double>> M(t->n_direction * 4);
+ std::vector<std::complex<double>> dM(t->n_direction * 16);
std::vector<double> dR(n_partial);
std::vector<double> dI(n_partial);
std::vector<unsigned int> dIndex(4 * n_partial); // 4 correlations
@@ -257,30 +260,30 @@ void grad_func(double *unknowns, double *grad, int m, void *adata) {
for (std::size_t dr = 0; dr < t->n_direction; ++dr) {
// Jones matrix for station P.
const double *Jp = &(unknowns[dr * t->n_station * 8 + p * 8]);
- const dcomplex Jp_00(Jp[0], Jp[1]);
- const dcomplex Jp_01(Jp[2], Jp[3]);
- const dcomplex Jp_10(Jp[4], Jp[5]);
- const dcomplex Jp_11(Jp[6], Jp[7]);
+ const std::complex<double> Jp_00(Jp[0], Jp[1]);
+ const std::complex<double> Jp_01(Jp[2], Jp[3]);
+ const std::complex<double> Jp_10(Jp[4], Jp[5]);
+ const std::complex<double> Jp_11(Jp[6], Jp[7]);
// Jones matrix for station Q, conjugated.
const double *Jq = &(unknowns[dr * t->n_station * 8 + q * 8]);
- const dcomplex Jq_00(Jq[0], -Jq[1]);
- const dcomplex Jq_01(Jq[2], -Jq[3]);
- const dcomplex Jq_10(Jq[4], -Jq[5]);
- const dcomplex Jq_11(Jq[6], -Jq[7]);
+ const std::complex<double> Jq_00(Jq[0], -Jq[1]);
+ const std::complex<double> Jq_01(Jq[2], -Jq[3]);
+ const std::complex<double> Jq_10(Jq[4], -Jq[5]);
+ const std::complex<double> Jq_11(Jq[6], -Jq[7]);
// Fetch model coherencies for the current direction.
- const dcomplex xx = t->model[dr][0];
- const dcomplex xy = t->model[dr][1];
- const dcomplex yx = t->model[dr][2];
- const dcomplex yy = t->model[dr][3];
+ const std::complex<double> xx = t->model[dr][0];
+ const std::complex<double> xy = t->model[dr][1];
+ const std::complex<double> yx = t->model[dr][2];
+ const std::complex<double> yy = t->model[dr][3];
// Precompute terms involving conj(Jq) and the model
// visibilities.
- const dcomplex Jq_00xx_01xy = Jq_00 * xx + Jq_01 * xy;
- const dcomplex Jq_00yx_01yy = Jq_00 * yx + Jq_01 * yy;
- const dcomplex Jq_10xx_11xy = Jq_10 * xx + Jq_11 * xy;
- const dcomplex Jq_10yx_11yy = Jq_10 * yx + Jq_11 * yy;
+ const std::complex<double> Jq_00xx_01xy = Jq_00 * xx + Jq_01 * xy;
+ const std::complex<double> Jq_00yx_01yy = Jq_00 * yx + Jq_01 * yy;
+ const std::complex<double> Jq_10xx_11xy = Jq_10 * xx + Jq_11 * xy;
+ const std::complex<double> Jq_10yx_11yy = Jq_10 * yx + Jq_11 * yy;
// Precompute (Jp x conj(Jq)) * vec(model), where 'x'
// denotes the Kronecker product. This is the model
@@ -322,75 +325,75 @@ void grad_func(double *unknowns, double *grad, int m, void *adata) {
if (!(t->flag[cr])) {
const double mwt = static_cast<double>(t->weight[cr]);
for (std::size_t tg = 0; tg < t->n_direction; ++tg) {
- dcomplex visibility(0.0, 0.0);
+ std::complex<double> visibility(0.0, 0.0);
for (std::size_t dr = 0; dr < t->n_direction; ++dr) {
// Look-up mixing weight.
- const dcomplex mix_weight = *(t->mix);
+ const std::complex<double> mix_weight = *(t->mix);
// Weight model visibility.
visibility += mix_weight * M[dr * 4 + cr];
// Compute weighted partial derivatives.
- dcomplex derivative(0.0, 0.0);
+ std::complex<double> derivative(0.0, 0.0);
derivative = mix_weight * dM[dr * 16 + cr * 4];
- dR[dr * 8] = real(derivative); // for cr==0: Re(d/dRe(p_00)))
- dI[dr * 8] = imag(derivative); // for cr==0: Re(d/dIm(p_00)))
+ dR[dr * 8] = derivative.real(); // for cr==0: Re(d/dRe(p_00)))
+ dI[dr * 8] = derivative.imag(); // for cr==0: Re(d/dIm(p_00)))
dR[dr * 8 + 1] =
- -imag(derivative); // for cr==0: Im(d/dRe(p_00)))
+ -derivative.imag(); // for cr==0: Im(d/dRe(p_00)))
dI[dr * 8 + 1] =
- real(derivative); // for cr==0: Im(d/dIm(p_00)))
+ derivative.real(); // for cr==0: Im(d/dIm(p_00)))
derivative = mix_weight * dM[dr * 16 + cr * 4 + 1];
dR[dr * 8 + 2] =
- real(derivative); // for cr==0: Re(d/dRe(p_01)))
+ derivative.real(); // for cr==0: Re(d/dRe(p_01)))
dI[dr * 8 + 2] =
- imag(derivative); // for cr==0: Re(d/dIm(p_01)))
+ derivative.imag(); // for cr==0: Re(d/dIm(p_01)))
dR[dr * 8 + 3] =
- -imag(derivative); // for cr==0: Im(d/dRe(p_01)))
+ -derivative.imag(); // for cr==0: Im(d/dRe(p_01)))
dI[dr * 8 + 3] =
- real(derivative); // for cr==0: Im(d/dIm(p_01)))
+ derivative.real(); // for cr==0: Im(d/dIm(p_01)))
derivative = mix_weight * dM[dr * 16 + cr * 4 + 2];
dR[dr * 8 + 4] =
- real(derivative); // for cr==0: Re(d/dRe(q_00)))
+ derivative.real(); // for cr==0: Re(d/dRe(q_00)))
dI[dr * 8 + 4] =
- imag(derivative); // for cr==0: Re(d/dIm(q_00)))
+ derivative.imag(); // for cr==0: Re(d/dIm(q_00)))
dR[dr * 8 + 5] =
- imag(derivative); // for cr==0: Im(d/dRe(q_00)))
+ derivative.imag(); // for cr==0: Im(d/dRe(q_00)))
dI[dr * 8 + 5] =
- -real(derivative); // for cr==0: Im(d/dIm(q_00)))
+ -derivative.real(); // for cr==0: Im(d/dIm(q_00)))
derivative = mix_weight * dM[dr * 16 + cr * 4 + 3];
dR[dr * 8 + 6] =
- real(derivative); // for cr==0: Re(d/dRe(q_01)))
+ derivative.real(); // for cr==0: Re(d/dRe(q_01)))
dI[dr * 8 + 6] =
- imag(derivative); // for cr==0: Re(d/dIm(q_01)))
+ derivative.imag(); // for cr==0: Re(d/dIm(q_01)))
dR[dr * 8 + 7] =
- imag(derivative); // for cr==0: Im(d/dRe(q_01)))
+ derivative.imag(); // for cr==0: Im(d/dRe(q_01)))
dI[dr * 8 + 7] =
- -real(derivative); // for cr==0: Im(d/dIm(q_01)))
+ -derivative.real(); // for cr==0: Im(d/dIm(q_01)))
// Move to next source direction.
t->mix.forward(1);
} // Source directions.
// Compute the residual.
- dcomplex residual =
- static_cast<dcomplex>(t->data[tg][cr]) - visibility;
+ const std::complex<double> residual =
+ std::complex<double>{t->data[tg][cr]} - visibility;
// accumulate gradient (for this correlation 'cr')
// For reference, gradient for Gaussian noise is:
// grad[dIndex[cr * n_partial + ci]] +=
- // 2.0 * dR[ci] * mwt * real(residual);
+ // 2.0 * dR[ci] * mwt * residual.real();
// grad[dIndex[cr * n_partial + ci]] +=
- // 2.0 * dI[ci] * mwt * imag(residual);
+ // 2.0 * dI[ci] * mwt * residual.imag();
for (std::size_t ci = 0; ci < n_partial; ci++) {
grad[dIndex[cr * n_partial + ci]] -=
- 2.0 * dR[ci] * mwt * real(residual) /
- (t->robust_nu + mwt * real(residual) * real(residual));
+ 2.0 * dR[ci] * mwt * residual.real() /
+ (t->robust_nu + mwt * residual.real() * residual.real());
grad[dIndex[cr * n_partial + ci]] -=
- 2.0 * dI[ci] * mwt * imag(residual) /
- (t->robust_nu + mwt * imag(residual) * imag(residual));
+ 2.0 * dI[ci] * mwt * residual.imag() /
+ (t->robust_nu + mwt * residual.imag() * residual.imag());
}
// Move to next target direction.
@@ -456,10 +459,10 @@ void grad_func(double *unknowns, double *grad, int m, void *adata) {
bool estimate(std::size_t n_direction, std::size_t n_station,
std::size_t n_baseline, std::size_t n_channel,
const_cursor<Baseline> baselines,
- std::vector<const_cursor<fcomplex>> data,
- std::vector<const_cursor<dcomplex>> model,
+ std::vector<const_cursor<std::complex<float>>> data,
+ std::vector<const_cursor<std::complex<double>>> model,
const_cursor<bool> flag, const_cursor<float> weight,
- const_cursor<dcomplex> mix, double *unknowns,
+ const_cursor<std::complex<double>> mix, double *unknowns,
std::size_t lbfgs_mem, double robust_nu, std::size_t max_iter) {
LBFGSData t(n_direction, n_station, n_baseline, n_channel, baselines, data,
model, flag, weight, mix, robust_nu);
diff --git a/base/Simulate.h b/base/Simulate.h
index a3c179b3c62e35f0f3ba76df1b1e2d7ab6be2302..05a89cc530ba4eb572881bb0338233c707b5678d 100644
--- a/base/Simulate.h
+++ b/base/Simulate.h
@@ -119,7 +119,7 @@ void simulate(const Direction& reference,
const std::shared_ptr<const Patch>& patch, size_t nStation,
size_t nBaseline, size_t nChannel,
const_cursor<Baseline> baselines, const_cursor<double> freq,
- const_cursor<double> uvw, cursor<dcomplex> buffer);
+ const_cursor<double> uvw, cursor<std::complex<double>> buffer);
} // namespace base
} // namespace dp3
diff --git a/base/SubtractMixed.cc b/base/SubtractMixed.cc
index ea40d398a530aee141f06a4e471b10a88ad0c659..3b1b2d5183c8b88c051d96d346566e6056de981d 100644
--- a/base/SubtractMixed.cc
+++ b/base/SubtractMixed.cc
@@ -12,8 +12,10 @@ namespace dp3 {
namespace base {
void subtract(size_t nBaseline, size_t nChannel,
- const_cursor<Baseline> baselines, cursor<fcomplex> data,
- const_cursor<dcomplex> model, const_cursor<dcomplex> weight) {
+ const_cursor<Baseline> baselines,
+ cursor<std::complex<float>> data,
+ const_cursor<std::complex<double>> model,
+ const_cursor<std::complex<double>> weight) {
for (size_t bl = 0; bl < nBaseline; ++bl) {
const size_t p = baselines->first;
const size_t q = baselines->second;
@@ -21,19 +23,19 @@ void subtract(size_t nBaseline, size_t nChannel,
if (p != q) {
for (size_t ch = 0; ch < nChannel; ++ch) {
// Subtract weighted model from data.
- *data -= static_cast<fcomplex>((*weight) * (*model));
+ *data -= std::complex<float>((*weight) * (*model));
++weight;
++model;
++data;
- *data -= static_cast<fcomplex>((*weight) * (*model));
+ *data -= std::complex<float>((*weight) * (*model));
++weight;
++model;
++data;
- *data -= static_cast<fcomplex>((*weight) * (*model));
+ *data -= std::complex<float>((*weight) * (*model));
++weight;
++model;
++data;
- *data -= static_cast<fcomplex>((*weight) * (*model));
+ *data -= std::complex<float>((*weight) * (*model));
++weight;
++model;
++data;
diff --git a/base/SubtractMixed.h b/base/SubtractMixed.h
index 362e682a6663919309a109b9192123a73031df18..6377540747e58ab7d3f4d5537ab56f2adedb0127 100644
--- a/base/SubtractMixed.h
+++ b/base/SubtractMixed.h
@@ -34,8 +34,10 @@ namespace base {
/// A cursor for a 3-D buffer of mixing weight of shape
/// (\p nBaseline, \p nChannel, 4).
void subtract(size_t nBaseline, size_t nChannel,
- const_cursor<Baseline> baselines, cursor<fcomplex> data,
- const_cursor<dcomplex> model, const_cursor<dcomplex> weight);
+ const_cursor<Baseline> baselines,
+ cursor<std::complex<float>> data,
+ const_cursor<std::complex<double>> model,
+ const_cursor<std::complex<double>> weight);
} // namespace base
} // namespace dp3
diff --git a/blob/BlobOStream.cc b/blob/BlobOStream.cc
index 000d34320f9bf936d410dec0f739509e7d54ef47..0744e6ece7d932dbe9e29a3d7f0761c385211935 100644
--- a/blob/BlobOStream.cc
+++ b/blob/BlobOStream.cc
@@ -137,11 +137,11 @@ BlobOStream& BlobOStream::operator<<(const double& var) {
putBuf(&var, sizeof(var));
return *this;
}
-BlobOStream& BlobOStream::operator<<(const fcomplex& var) {
+BlobOStream& BlobOStream::operator<<(const std::complex<float>& var) {
putBuf(&var, sizeof(var));
return *this;
}
-BlobOStream& BlobOStream::operator<<(const dcomplex& var) {
+BlobOStream& BlobOStream::operator<<(const std::complex<double>& var) {
putBuf(&var, sizeof(var));
return *this;
}
@@ -201,11 +201,11 @@ void BlobOStream::put(const float* values, uint64_t nrval) {
void BlobOStream::put(const double* values, uint64_t nrval) {
putBuf(values, nrval * sizeof(double));
}
-void BlobOStream::put(const fcomplex* values, uint64_t nrval) {
- putBuf(values, nrval * sizeof(fcomplex));
+void BlobOStream::put(const std::complex<float>* values, uint64_t nrval) {
+ putBuf(values, nrval * sizeof(std::complex<float>));
}
-void BlobOStream::put(const dcomplex* values, uint64_t nrval) {
- putBuf(values, nrval * sizeof(dcomplex));
+void BlobOStream::put(const std::complex<double>* values, uint64_t nrval) {
+ putBuf(values, nrval * sizeof(std::complex<double>));
}
void BlobOStream::put(const std::string* values, uint64_t nrval) {
for (uint64_t i = 0; i < nrval; i++) {
diff --git a/blob/BlobOStream.h b/blob/BlobOStream.h
index 3264ec9969d299f82c6199d5b1e6300dc4ab9eb3..b6ec57c2d20a0d473e7aa647010af2bf5b1414b1 100644
--- a/blob/BlobOStream.h
+++ b/blob/BlobOStream.h
@@ -14,9 +14,6 @@
#include <complex>
#include <cstring>
-typedef std::complex<double> dcomplex;
-typedef std::complex<float> fcomplex;
-
namespace dp3 {
namespace blob {
@@ -83,12 +80,6 @@ class BlobOStream {
/// Put a single value.
/// A string will be stored as a length followed by the characters.
- /// Note that a function is defined for the standard complex class
- /// and for the types fcomplex, dcomplex, i16complex and u16complex.
- /// This should be fine for the case where fcomplex is the builtin
- /// complex type and the case where it is the standard class.
- /// In the first case the fcomplex function is a separate function,
- /// in the second case a specialisation of the templated function.
/// @{
BlobOStream& operator<<(const bool& value);
BlobOStream& operator<<(const char& value);
@@ -102,13 +93,8 @@ class BlobOStream {
BlobOStream& operator<<(const uint64_t& value);
BlobOStream& operator<<(const float& value);
BlobOStream& operator<<(const double& value);
- template <class T>
- BlobOStream& operator<<(const std::complex<T>& value);
- // BlobOStream& operator<< (const i4complex& value);
- // BlobOStream& operator<< (const i16complex& value);
- // BlobOStream& operator<< (const u16complex& value);
- BlobOStream& operator<<(const fcomplex& value);
- BlobOStream& operator<<(const dcomplex& value);
+ BlobOStream& operator<<(const std::complex<float>& value);
+ BlobOStream& operator<<(const std::complex<double>& value);
BlobOStream& operator<<(const std::string& value);
BlobOStream& operator<<(const char* value);
/// @}
@@ -128,13 +114,8 @@ class BlobOStream {
void put(const uint64_t* values, uint64_t nrval);
void put(const float* values, uint64_t nrval);
void put(const double* values, uint64_t nrval);
- template <class T>
- void put(const std::complex<T>* values, uint64_t nrval);
- // void put (const i4complex* values, uint64_t nrval);
- // void put (const i16complex* values, uint64_t nrval);
- // void put (const u16complex* values, uint64_t nrval);
- void put(const fcomplex* values, uint64_t nrval);
- void put(const dcomplex* values, uint64_t nrval);
+ void put(const std::complex<float>* values, uint64_t nrval);
+ void put(const std::complex<double>* values, uint64_t nrval);
void put(const std::string* values, uint64_t nrval);
/// @}
@@ -217,16 +198,6 @@ inline unsigned int BlobOStream::putStart(const char* objectType,
inline int64_t BlobOStream::tellPos() const { return itsStream->tellPos(); }
-template <class T>
-inline BlobOStream& BlobOStream::operator<<(const std::complex<T>& value) {
- putBuf(&value, sizeof(value));
- return *this;
-}
-template <class T>
-inline void BlobOStream::put(const std::complex<T>* values, uint64_t nrval) {
- putBuf(values, nrval * sizeof(std::complex<T>));
-}
-
template <class T>
inline void BlobOStream::put(const std::vector<T>& vec) {
operator<<(uint64_t(vec.size()));
diff --git a/common/TypeNames.cc b/common/TypeNames.cc
index b617b8d205f37dde5f1e833127afd81ecc4df0eb..6f79a8e015ea23f0f309f03cfc31bdd3bb7a67eb 100644
--- a/common/TypeNames.cc
+++ b/common/TypeNames.cc
@@ -80,6 +80,7 @@ const std::string& typeName(const std::complex<float>*) {
static std::string str("fcomplex");
return str;
}
+
const std::string& typeName(const std::complex<double>*) {
static std::string str("dcomplex");
return str;
diff --git a/steps/Demixer.cc b/steps/Demixer.cc
index 503227f7b7c69d2e9b19c263c30a10fce430c464..321a0de70981a2488129791180ba7bb3a85c77c5 100644
--- a/steps/Demixer.cc
+++ b/steps/Demixer.cc
@@ -870,8 +870,8 @@ namespace {
struct ThreadPrivateStorage {
std::vector<double> unknowns;
xt::xtensor<double, 2> uvw;
- std::vector<casacore::Cube<dcomplex>> model;
- casacore::Cube<dcomplex> model_subtr;
+ std::vector<casacore::Cube<std::complex<double>>> model;
+ casacore::Cube<std::complex<double>> model_subtr;
size_t count_converged;
};
@@ -977,17 +977,18 @@ void Demixer::demix() {
itsFactors[ts].shape(0));
const casacore::Array<casacore::DComplex> demix_factor(
shape, itsFactors[ts].data(), casacore::SHARE);
- base::const_cursor<dcomplex> cr_mix = base::casa_const_cursor(demix_factor);
+ base::const_cursor<std::complex<double>> cr_mix =
+ base::casa_const_cursor(demix_factor);
/// cout << "demixfactor "<<ts<<" = "<<itsFactors[ts]<<'\n';
- std::vector<base::const_cursor<fcomplex>> cr_data(nDr);
- std::vector<base::const_cursor<dcomplex>> cr_model(nDr);
+ std::vector<base::const_cursor<std::complex<float>>> cr_data(nDr);
+ std::vector<base::const_cursor<std::complex<double>>> cr_model(nDr);
for (size_t dr = 0; dr < nDr; ++dr) {
cr_data[dr] = base::casa_const_cursor(casacore::Cube<casacore::Complex>(
cube_shape, itsAvgResults[dr]->get()[ts]->GetData().data(),
casacore::SHARE));
- cr_model[dr] = base::const_cursor<dcomplex>(storage.model[dr].data(), 3,
- stride_model);
+ cr_model[dr] = base::const_cursor<std::complex<double>>(
+ storage.model[dr].data(), 3, stride_model);
}
const bool converged =
@@ -1021,8 +1022,8 @@ void Demixer::demix() {
ts_subtr != ts_subtr_end; ++ts_subtr) {
for (size_t dr = 0; dr < nDrSubtr; ++dr) {
// Re-use simulation used for estimating Jones matrices if possible.
- base::cursor<dcomplex> cr_model_subtr(storage.model[dr].data(), 3,
- stride_model);
+ base::cursor<std::complex<double>> cr_model_subtr(
+ storage.model[dr].data(), 3, stride_model);
// Re-simulate if required.
if (multiplier != 1 || nCh != nChSubtr) {
@@ -1051,12 +1052,12 @@ void Demixer::demix() {
storage.uvw.data());
// Zero the visibility buffer.
- storage.model_subtr = dcomplex();
+ storage.model_subtr = std::complex<double>();
// Simulate visibilities at the resolution of the residual.
size_t stride_model_subtr[3] = {1, nCr, nCr * nChSubtr};
- cr_model_subtr = base::cursor<dcomplex>(storage.model_subtr.data(), 3,
- stride_model_subtr);
+ cr_model_subtr = base::cursor<std::complex<double>>(
+ storage.model_subtr.data(), 3, stride_model_subtr);
base::Simulator simulator(itsPatchList[dr]->direction(), nSt,
itsBaselines,
@@ -1083,7 +1084,8 @@ void Demixer::demix() {
casacore::Cube<casacore::Complex> cube_residual(
cube_shape, itsAvgResultSubtr->get()[ts_subtr]->GetData().data(),
casacore::SHARE);
- base::cursor<fcomplex> cr_residual = base::casa_cursor(cube_residual);
+ base::cursor<std::complex<float>> cr_residual =
+ base::casa_cursor(cube_residual);
// Construct a cursor to iterate over a slice of the mixing matrix
// at the resolution of the residual. The "to" and "from" direction
@@ -1109,9 +1111,8 @@ void Demixer::demix() {
itsFactorsSubtr[ts_subtr].data(), casacore::SHARE);
IPosition offset(5, itsNDir - 1, dr, 0, 0, 0);
- base::const_cursor<dcomplex> cr_mix_subtr =
- base::const_cursor<dcomplex>(&(demix_factors_subtr(offset)), 3,
- stride_mix_subtr_slice);
+ base::const_cursor<std::complex<double>> cr_mix_subtr(
+ &(demix_factors_subtr(offset)), 3, stride_mix_subtr_slice);
// Subtract the source.
subtract(nBl, nChSubtr, cr_baseline, cr_residual, cr_model_subtr,
diff --git a/steps/OnePredict.cc b/steps/OnePredict.cc
index 147fa6d601807a278cb098072d4c45c0ec5c9d31..306a54c6c0c138f61ad7180c81e55d94efdc1922 100644
--- a/steps/OnePredict.cc
+++ b/steps/OnePredict.cc
@@ -55,7 +55,6 @@
#include <boost/algorithm/string/case_conv.hpp>
-using casacore::Cube;
using casacore::MDirection;
using casacore::MEpoch;
using casacore::MVDirection;
@@ -479,7 +478,8 @@ bool OnePredict::process(std::unique_ptr<DPBuffer> buffer) {
const casacore::IPosition shape(3, thread_buffer.shape(2),
thread_buffer.shape(1),
thread_buffer.shape(0));
- Cube<dcomplex> simulatedest(shape, thread_buffer.data(), casacore::SHARE);
+ casacore::Cube<std::complex<double>> simulatedest(
+ shape, thread_buffer.data(), casacore::SHARE);
simulators.emplace_back(phase_ref_, nSt, baselines_split[thread_index],
casacore::Vector<double>(info().chanFreqs()),
@@ -507,7 +507,8 @@ bool OnePredict::process(std::unique_ptr<DPBuffer> buffer) {
// Always use model.shape(), since it's equal to the patch model shape.
const casacore::IPosition shape(3, model.shape(2), model.shape(1),
model.shape(0));
- Cube<dcomplex> simulatedest(shape, model_data, casacore::SHARE);
+ casacore::Cube<std::complex<double>> simulatedest(shape, model_data,
+ casacore::SHARE);
simulators.emplace_back(phase_ref_, nSt, baselines_,
casacore::Vector<double>(info().chanFreqs()),