diff --git a/.gitlab-ci.common.yml b/.gitlab-ci.common.yml
index de3ffcc0b17b568820ddb1f3a886068457d2fa7e..45def8e44f4727400b64422308a26d794307feda 100644
--- a/.gitlab-ci.common.yml
+++ b/.gitlab-ci.common.yml
@@ -11,6 +11,9 @@ workflow:
when: never
- when: always
+variables:
+ GIT_SUBMODULE_STRATEGY: recursive
+
stages:
- versioning
- prepare
diff --git a/.gitmodules b/.gitmodules
index d34218b42bd3d311fb017efb5be9083961cb0b7c..58a5a02c14f3370ff3663470518ec6de97c7371a 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,3 +1,6 @@
[submodule "external/schaap-packaging"]
path = external/schaap-packaging
url = https://git.astron.nl/RD/schaap-packaging.git
+[submodule "external/aocommon"]
+ path = external/aocommon
+ url = https://gitlab.com/aroffringa/aocommon.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a176621e1310c3d1216a5bbdbf54df5857077de4..a8e6fed7edb1d51f71f74e2034f7f193178c09ec 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -30,6 +30,11 @@ set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED YES)
set(CMAKE_CXX_EXTENSIONS NO)
+include_directories(${CMAKE_SOURCE_DIR}/external/aocommon/include/)
+include(external/aocommon/CMake/FetchXTensor.cmake)
+include_directories(SYSTEM "${xtensor_SOURCE_DIR}/include")
+include_directories(SYSTEM "${xtl_SOURCE_DIR}/include")
+
# build shared or static libraries (default: shared)
if(BUILD_STATIC_LIBS)
set(BUILD_SHARED_LIBS FALSE)
diff --git a/external/aocommon b/external/aocommon
new file mode 160000
index 0000000000000000000000000000000000000000..c9b57adbb537a94a5dd5501697a31fa6bbb26aee
--- /dev/null
+++ b/external/aocommon
@@ -0,0 +1 @@
+Subproject commit c9b57adbb537a94a5dd5501697a31fa6bbb26aee
diff --git a/idg-api/BufferSet.cpp b/idg-api/BufferSet.cpp
index 9c7b4a78e4635d530a83f0bd231f48ea44c5001f..da15dd8a80e6ca11242ed8bc348d2ebdc78fe9e8 100644
--- a/idg-api/BufferSet.cpp
+++ b/idg-api/BufferSet.cpp
@@ -76,15 +76,13 @@ BufferSetImpl::BufferSetImpl(Type architecture)
m_nr_correlations(4),
m_nr_polarizations(4),
m_proxy(create_proxy(architecture)),
- m_shift(2),
+ m_shift({0, 0}),
m_get_image_watch(Stopwatch::create()),
m_set_image_watch(Stopwatch::create()),
m_avg_beam_watch(Stopwatch::create()),
m_plan_watch(Stopwatch::create()),
m_gridding_watch(Stopwatch::create()),
- m_degridding_watch(Stopwatch::create()) {
- m_shift.zero();
-}
+ m_degridding_watch(Stopwatch::create()) {}
BufferSetImpl::~BufferSetImpl() {
// Free all objects allocated via the proxy before destroying the proxy.
@@ -258,8 +256,8 @@ void BufferSetImpl::init(size_t size, float cell_size, float max_w,
std::cout << "nr_w_layers: " << m_nr_w_layers << std::endl;
#endif
- m_shift(0) = shiftl;
- m_shift(1) = shiftm;
+ m_shift[0] = shiftl;
+ m_shift[1] = shiftm;
m_kernel_size = taper_kernel_size + w_kernel_size + a_term_kernel_size;
diff --git a/idg-api/BufferSetImpl.h b/idg-api/BufferSetImpl.h
index eee20a38ce7bc98d8a7a53fcde46650b592d52fe..44190a4d04e1144d33d99c0fb5caff76607c4fa1 100644
--- a/idg-api/BufferSetImpl.h
+++ b/idg-api/BufferSetImpl.h
@@ -78,7 +78,7 @@ class BufferSetImpl : public virtual BufferSet {
float get_cell_size() const { return m_cell_size; }
float get_w_step() const { return m_w_step; }
- const idg::Array1D<float>& get_shift() const { return m_shift; }
+ const std::array<float, 2>& get_shift() const { return m_shift; }
float get_kernel_size() const { return m_kernel_size; }
const Array2D<float>& get_spheroidal() const { return m_spheroidal; }
@@ -123,7 +123,7 @@ class BufferSetImpl : public virtual BufferSet {
float m_cell_size;
float m_w_step;
int m_nr_w_layers;
- Array1D<float> m_shift;
+ std::array<float, 2> m_shift;
size_t m_size;
size_t m_padded_size;
float m_kernel_size;
diff --git a/idg-api/BulkDegridder.cpp b/idg-api/BulkDegridder.cpp
index 1434ccf316576b3927a4cf8bfc93c5a8c34be93f..f37d061c73332ab45d615761524557849481796f 100644
--- a/idg-api/BulkDegridder.cpp
+++ b/idg-api/BulkDegridder.cpp
@@ -58,19 +58,22 @@ void BulkDegridderImpl::compute_visibilities(
static const double kDefaultUVWFactors[3] = {1.0, 1.0, 1.0};
if (!uvw_factors) uvw_factors = kDefaultUVWFactors;
- auto bufferStationPairs =
- proxy.allocate_array1d<std::pair<unsigned int, unsigned int>>(
- nr_baselines);
+ auto bufferStationPairs_tensor =
+ proxy.allocate_tensor<std::pair<unsigned int, unsigned int>, 1>(
+ {nr_baselines});
+ auto bufferStationPairs = bufferStationPairs_tensor.Span();
+ bufferStationPairs.fill(
+ std::pair<unsigned int, unsigned int>(nr_stations_, nr_stations_));
std::vector<int> baseline_map;
- auto bufferUVW =
- proxy.allocate_array2d<UVW<float>>(nr_baselines, nr_timesteps);
+ auto bufferUVW_tensor =
+ proxy.allocate_tensor<UVW<float>, 2>({nr_baselines, nr_timesteps});
+ auto bufferUVW = bufferUVW_tensor.Span();
auto bufferVisibilities = proxy.allocate_array4d<std::complex<float>>(
nr_baselines, nr_timesteps, frequencies_.size(), nr_correlations);
// The proxy does not touch visibilities for out-of-bound uvw coordinates.
// Since we copy all visibilities to the caller, initialize them to zero.
bufferVisibilities.zero();
- bufferStationPairs.init({nr_stations_, nr_stations_});
baseline_map.reserve(nr_baselines);
for (size_t bl = 0; bl < nr_baselines; ++bl) {
const size_t ant1 = antennas1[bl];
@@ -96,8 +99,10 @@ void BulkDegridderImpl::compute_visibilities(
// The proxy expects the number of time steps as last value.
std::vector<unsigned int> local_aterm_offsets = aterm_offsets;
local_aterm_offsets.push_back(nr_timesteps);
- const Array1D<unsigned int> aterm_offsets_array(local_aterm_offsets.data(),
- local_aterm_offsets.size());
+ const std::array<size_t, 1> local_aterm_offsets_shape{
+ local_aterm_offsets.size()};
+ auto aterm_offsets_span = aocommon::xt::CreateSpan(local_aterm_offsets.data(),
+ local_aterm_offsets_shape);
// If aterms is empty, create default values and update the pointer.
std::vector<Matrix2x2<std::complex<float>>> default_aterms;
@@ -113,8 +118,8 @@ void BulkDegridderImpl::compute_visibilities(
auto* local_aterms = reinterpret_cast<Matrix2x2<std::complex<float>>*>(
const_cast<std::complex<float>*>(aterms));
const Array4D<Matrix2x2<std::complex<float>>> aterms_array(
- local_aterms, aterm_offsets_array.get_x_dim() - 1, nr_stations_,
- subgridsize, subgridsize);
+ local_aterms, aterm_offsets_span.size() - 1, nr_stations_, subgridsize,
+ subgridsize);
// Set Plan options
Plan::Options options;
@@ -126,15 +131,18 @@ void BulkDegridderImpl::compute_visibilities(
// Create plan
bufferset_.get_watch(BufferSetImpl::Watch::kPlan).Start();
+ const std::array<size_t, 1> frequencies_shape{frequencies_.size()};
+ auto frequencies_span =
+ aocommon::xt::CreateSpan(frequencies_.data(), frequencies_shape);
std::unique_ptr<Plan> plan =
- proxy.make_plan(bufferset_.get_kernel_size(), frequencies_, bufferUVW,
- bufferStationPairs, aterm_offsets_array, options);
+ proxy.make_plan(bufferset_.get_kernel_size(), frequencies_span, bufferUVW,
+ bufferStationPairs, aterm_offsets_span, options);
bufferset_.get_watch(BufferSetImpl::Watch::kPlan).Pause();
// Run degridding
bufferset_.get_watch(BufferSetImpl::Watch::kDegridding).Start();
proxy.degridding(*plan, frequencies_, bufferVisibilities, bufferUVW,
- bufferStationPairs, aterms_array, aterm_offsets_array,
+ bufferStationPairs, aterms_array, aterm_offsets_span,
bufferset_.get_spheroidal());
bufferset_.get_watch(BufferSetImpl::Watch::kDegridding).Pause();
diff --git a/idg-api/DegridderBuffer.cpp b/idg-api/DegridderBuffer.cpp
index 8e09f261923c1f9920cdbd332a492524bb7d12f0..6c4605bfbdc081edafb364d64f9b2529d9f90bdd 100644
--- a/idg-api/DegridderBuffer.cpp
+++ b/idg-api/DegridderBuffer.cpp
@@ -103,6 +103,10 @@ void DegridderBufferImpl::flush() {
m_aterms.data(), m_aterm_offsets_array.get_x_dim() - 1, m_nrStations,
subgridsize, subgridsize);
+ const std::array<size_t, 1> aterm_offsets_shape{m_aterm_offsets_array.size()};
+ auto aterm_offsets_span = aocommon::xt::CreateSpan(
+ m_aterm_offsets_array.data(), aterm_offsets_shape);
+
proxy::Proxy& proxy = m_bufferset.get_proxy();
// Set Plan options
@@ -115,9 +119,18 @@ void DegridderBufferImpl::flush() {
// Create plan
m_bufferset.get_watch(BufferSetImpl::Watch::kPlan).Start();
+ const std::array<size_t, 1> frequencies_shape{m_frequencies.size()};
+ auto frequencies_span =
+ aocommon::xt::CreateSpan(m_frequencies.data(), frequencies_shape);
+ const std::array<size_t, 1> station_pairs_shape{m_bufferStationPairs.size()};
+ const std::array<size_t, 2> uvw_shape{m_bufferUVW.get_y_dim(),
+ m_bufferUVW.get_x_dim()};
+ auto uvw_span = aocommon::xt::CreateSpan(m_bufferUVW.data(), uvw_shape);
+ auto station_pairs_span = aocommon::xt::CreateSpan(
+ m_bufferStationPairs.data(), station_pairs_shape);
std::unique_ptr<Plan> plan =
- proxy.make_plan(m_bufferset.get_kernel_size(), m_frequencies, m_bufferUVW,
- m_bufferStationPairs, m_aterm_offsets_array, options);
+ proxy.make_plan(m_bufferset.get_kernel_size(), frequencies_span, uvw_span,
+ station_pairs_span, aterm_offsets_span, options);
m_bufferset.get_watch(BufferSetImpl::Watch::kPlan).Pause();
// Run degridding
diff --git a/idg-api/GridderBuffer.cpp b/idg-api/GridderBuffer.cpp
index 084c7c33f83f4383c193e32171c18ac747d50645..b5f7a4cf3f237866c77fdcb3712fa02cdb1bd69f 100644
--- a/idg-api/GridderBuffer.cpp
+++ b/idg-api/GridderBuffer.cpp
@@ -139,6 +139,10 @@ void GridderBufferImpl::flush_thread_worker() {
m_nrStations, subgridsize, subgridsize);
}
+ const std::array<size_t, 1> aterm_offsets_shape{m_aterm_offsets_array.size()};
+ auto aterm_offsets_span = aocommon::xt::CreateSpan(
+ m_aterm_offsets_array.data(), aterm_offsets_shape);
+
proxy::Proxy& proxy = m_bufferset.get_proxy();
// Set Plan options
@@ -151,9 +155,18 @@ void GridderBufferImpl::flush_thread_worker() {
// Create plan
m_bufferset.get_watch(BufferSetImpl::Watch::kPlan).Start();
- std::unique_ptr<Plan> plan = proxy.make_plan(
- m_bufferset.get_kernel_size(), m_frequencies, m_bufferUVW2,
- m_bufferStationPairs2, m_aterm_offsets_array, options);
+ const std::array<size_t, 1> frequencies_shape{m_frequencies.size()};
+ auto frequencies_span =
+ aocommon::xt::CreateSpan(m_frequencies.data(), frequencies_shape);
+ const std::array<size_t, 1> station_pairs_shape{m_bufferStationPairs.size()};
+ const std::array<size_t, 2> uvw_shape{m_bufferUVW2.get_y_dim(),
+ m_bufferUVW2.get_x_dim()};
+ auto uvw_span = aocommon::xt::CreateSpan(m_bufferUVW2.data(), uvw_shape);
+ auto station_pairs_span = aocommon::xt::CreateSpan(
+ m_bufferStationPairs2.data(), station_pairs_shape);
+ std::unique_ptr<Plan> plan =
+ proxy.make_plan(m_bufferset.get_kernel_size(), frequencies_span, uvw_span,
+ station_pairs_span, aterm_offsets_span, options);
m_bufferset.get_watch(BufferSetImpl::Watch::kPlan).Pause();
// Run gridding
diff --git a/idg-bin/CMakeLists.txt b/idg-bin/CMakeLists.txt
index 02cd98c066740ec34deacbdaf7c63017c4199ddf..a6f6b380f9a84c46c356c051a23a7706a5334c15 100644
--- a/idg-bin/CMakeLists.txt
+++ b/idg-bin/CMakeLists.txt
@@ -10,7 +10,6 @@ list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/modules)
# options
option(BUILD_LIB_CPU "Build CPU libraries" ON)
option(BUILD_LIB_CUDA "Build CUDA libraries" OFF)
-option(USE_DUMMY_VISIBILITIES "Run C++ examples using dummy visibilities" ON)
# Compiler settings:
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
diff --git a/idg-bin/examples/cxx/common/common.h b/idg-bin/examples/cxx/common/common.h
index 45d06f8d08b68fe723d32955f759ccd1df75057e..c7930899348d3294be8d4f21184e523ca070846a 100644
--- a/idg-bin/examples/cxx/common/common.h
+++ b/idg-bin/examples/cxx/common/common.h
@@ -15,14 +15,6 @@
#include "idg-cpu.h"
#include "idg-util.h" // Data init routines
-/*
- * Visibilities initializion
- * 0, use simulated visibilities for some point sources near the centre of the
- * image 1, use fixed (dummy) visibilities to speed-up initialization of the
- * data
- */
-#define USE_DUMMY_VISIBILITIES 1
-
using namespace std;
std::tuple<int, int, int, int, int, int, int, int, int, int, int, float, bool,
@@ -240,22 +232,22 @@ void run() {
// Allocate and initialize static data structures
clog << ">>> Initialize data structures" << endl;
-#if USE_DUMMY_VISIBILITIES
- idg::Array4D<std::complex<float>> visibilities_ = idg::get_dummy_visibilities(
- proxy, nr_baselines, nr_timesteps, nr_channels, nr_correlations);
-#endif
- idg::Array4D<idg::Matrix2x2<std::complex<float>>> aterms =
+ aocommon::xt::Span<std::complex<float>, 4> visibilities =
+ idg::get_dummy_visibilities(proxy, nr_baselines, nr_timesteps,
+ nr_channels, nr_correlations);
+ aocommon::xt::Span<idg::Matrix2x2<std::complex<float>>, 4> aterms =
idg::get_identity_aterms(proxy, nr_timeslots, nr_stations, subgrid_size,
subgrid_size);
- idg::Array1D<unsigned int> aterm_offsets =
+ aocommon::xt::Span<unsigned int, 1> aterm_offsets =
idg::get_example_aterm_offsets(proxy, nr_timeslots, nr_timesteps);
- idg::Array2D<float> spheroidal =
+ aocommon::xt::Span<float, 2> spheroidal =
idg::get_example_spheroidal(proxy, subgrid_size, subgrid_size);
+
auto grid =
proxy.allocate_grid(nr_w_layers, nr_polarizations, grid_size, grid_size);
grid->zero();
- idg::Array1D<float> shift = idg::get_zero_shift();
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines =
+ std::array<float, 2> shift{0.0f, 0.0f};
+ aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1> baselines =
idg::get_example_baselines(proxy, nr_stations, nr_baselines);
clog << endl;
@@ -272,13 +264,9 @@ void run() {
bool disable_degridding = getenv("DISABLE_DEGRIDDING");
bool disable_fft = getenv("DISABLE_FFT");
- // Spectral line imaging
- bool simulate_spectral_line = getenv("SPECTRAL_LINE");
-
// Plan options
idg::Plan::Options options;
options.plan_strict = true;
- options.simulate_spectral_line = simulate_spectral_line;
options.max_nr_timesteps_per_subgrid = 128;
options.max_nr_channels_per_subgrid = 8;
options.mode = stokes_i_only ? idg::Plan::Mode::STOKES_I_ONLY
@@ -304,13 +292,16 @@ void run() {
std::ceil((float)total_nr_timesteps / nr_timesteps);
for (unsigned int t = 0; t < nr_time_blocks; t++) {
int time_offset = t * nr_timesteps;
- int current_nr_timesteps = total_nr_timesteps - time_offset < nr_timesteps
- ? total_nr_timesteps - time_offset
- : nr_timesteps;
+ const size_t current_nr_timesteps =
+ total_nr_timesteps - time_offset < nr_timesteps
+ ? total_nr_timesteps - time_offset
+ : nr_timesteps;
// Initalize UVW coordiantes
- auto uvw = proxy.allocate_array2d<idg::UVW<float>>(nr_baselines,
- current_nr_timesteps);
+ aocommon::xt::Span<idg::UVW<float>, 2> uvw =
+ proxy.allocate_span<idg::UVW<float>, 2>(
+ {nr_baselines, current_nr_timesteps});
+
data.get_uvw(uvw, 0, time_offset, integration_time);
// Iterate all channel blocks
@@ -325,22 +316,10 @@ void run() {
clog << ">>>" << endl;
// Initialize frequency data
- idg::Array1D<float> frequencies =
- proxy.allocate_array1d<float>(nr_channels);
+ aocommon::xt::Span<float, 1> frequencies =
+ proxy.allocate_span<float, 1>({nr_channels});
data.get_frequencies(frequencies, image_size, channel_offset);
-// Initialize visibilities
-#if !USE_DUMMY_VISIBILITIES
- idg::Array4D<std::complex<float>> visibilities_ =
- idg::get_example_visibilities(proxy, uvw, frequencies, image_size,
- grid_size, nr_correlations);
-#endif
-
- auto nr_channels_ = simulate_spectral_line ? 1 : nr_channels;
- idg::Array4D<std::complex<float>> visibilities(
- visibilities_.data(), nr_baselines, current_nr_timesteps,
- nr_channels_, nr_correlations);
-
// Create plan
if (plans.size() == 0 || cycle == 0) {
plans.emplace_back(proxy.make_plan(kernel_size, frequencies, uvw,
diff --git a/idg-bin/examples/cxx/distributed/common/common.h b/idg-bin/examples/cxx/distributed/common/common.h
index adf720acc2d6c320fe95e03dea407ae0f326e0d3..8df1c479f56c4c57dacca9107e594b913fdb56e4 100644
--- a/idg-bin/examples/cxx/distributed/common/common.h
+++ b/idg-bin/examples/cxx/distributed/common/common.h
@@ -147,14 +147,15 @@ float receive_float(int src = 0) {
}
template <typename T>
-void send_array(int dst, T& array) {
- MPI_Send(array.data(), array.bytes(), MPI_BYTE, dst, 0, MPI_COMM_WORLD);
+void send_tensor(int dst, T& tensor) {
+ MPI_Send(tensor.data(), tensor.size() * sizeof(T), MPI_BYTE, dst, 0,
+ MPI_COMM_WORLD);
}
template <typename T>
-void receive_array(int src, T& array) {
- MPI_Recv(array.data(), array.bytes(), MPI_BYTE, src, 0, MPI_COMM_WORLD,
- MPI_STATUS_IGNORE);
+void receive_tensor(int src, T& tensor) {
+ MPI_Recv(tensor.data(), tensor.size() * sizeof(T), MPI_BYTE, src, 0,
+ MPI_COMM_WORLD, MPI_STATUS_IGNORE);
}
void send_bytes(int dst, void* buf, size_t bytes) {
@@ -249,7 +250,7 @@ void reduce_grids(std::shared_ptr<idg::Grid> grid, unsigned int rank,
unsigned int nr_polarizations = grid->get_z_dim();
unsigned int grid_size = grid->get_y_dim();
- idg::Array2D<std::complex<float>> tmp(grid_size, grid_size);
+ xt::xtensor<std::complex<float>, 2> tmp({grid_size, grid_size});
size_t sizeof_row = grid_size * sizeof(std::complex<float>);
#pragma omp parallel
@@ -260,7 +261,7 @@ void reduce_grids(std::shared_ptr<idg::Grid> grid, unsigned int rank,
if (omp_get_thread_num() == 0) {
MPIRequestList requests;
for (unsigned int y = 0; y < grid_size; y++) {
- requests.create()->receive(tmp.data(y, 0), sizeof_row, i + rank);
+ requests.create()->receive(&tmp(y, 0), sizeof_row, i + rank);
}
requests.wait();
}
@@ -271,12 +272,12 @@ void reduce_grids(std::shared_ptr<idg::Grid> grid, unsigned int rank,
#pragma omp for
for (unsigned int y = 0; y < grid_size; y++)
for (unsigned int x = 0; x < grid_size; x++) {
- grid_(w, pol, y, x) += *tmp.data(y, x);
+ grid_(w, pol, y, x) += tmp(y, x);
}
} else if (rank < (2 * i) && omp_get_thread_num() == 0) {
MPIRequestList requests;
for (unsigned int y = 0; y < grid_size; y++) {
- requests.create()->send(tmp.data(y, 0), sizeof_row, rank - i);
+ requests.create()->send(&tmp(y, 0), sizeof_row, rank - i);
}
}
} // end for i
@@ -383,8 +384,12 @@ void run_master() {
send_string(dst, layout_file);
}
+ // Initialize proxy
+ ProxyType proxy;
+
// Initialize frequency data for master
- idg::Array1D<float> frequencies(nr_channels);
+ aocommon::xt::Span<float, 1> frequencies =
+ proxy.allocate_span<float, 1>({nr_channels});
data.get_frequencies(frequencies, image_size);
// Distribute frequencies to workers
@@ -392,10 +397,11 @@ void run_master() {
// take a channel offset into account when initializing
// frequencies.
for (int dst = 1; dst < world_size; dst++) {
- int channel_offset = (dst - 1) * nr_channels;
- idg::Array1D<float> frequencies_(nr_channels);
+ const int channel_offset = (dst - 1) * nr_channels;
+ aocommon::xt::Span<float, 1> frequencies_ =
+ proxy.allocate_span<float, 1>({nr_channels});
data.get_frequencies(frequencies_, image_size, channel_offset);
- send_array(dst, frequencies_);
+ send_tensor(dst, frequencies_);
}
// Distribute data
@@ -403,21 +409,18 @@ void run_master() {
send_bytes(dst, &data, sizeof(data));
}
- // Initialize proxy
- ProxyType proxy;
-
// Allocate and initialize static data structures
- idg::Array4D<idg::Matrix2x2<std::complex<float>>> aterms =
+ aocommon::xt::Span<idg::Matrix2x2<std::complex<float>>, 4> aterms =
idg::get_identity_aterms(proxy, nr_timeslots, nr_stations, subgrid_size,
subgrid_size);
- idg::Array1D<unsigned int> aterm_offsets =
+ aocommon::xt::Span<unsigned int, 1> aterm_offsets =
idg::get_example_aterm_offsets(proxy, nr_timeslots, nr_timesteps);
- idg::Array2D<float> spheroidal =
+ aocommon::xt::Span<float, 2> spheroidal =
idg::get_example_spheroidal(proxy, subgrid_size, subgrid_size);
- auto grid =
+ std::shared_ptr<idg::Grid> grid =
proxy.allocate_grid(nr_w_layers, nr_correlations, grid_size, grid_size);
- idg::Array1D<float> shift = idg::get_zero_shift();
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines =
+ std::array<float, 2> shift{0.0f, 0.0f};
+ aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1> baselines =
idg::get_example_baselines(proxy, nr_stations, nr_baselines);
// Plan options
@@ -427,10 +430,12 @@ void run_master() {
// Buffers for input data
unsigned int nr_time_blocks =
std::ceil((float)total_nr_timesteps / nr_timesteps);
- idg::Array3D<idg::UVW<float>> uvws = proxy.allocate_array3d<idg::UVW<float>>(
- nr_time_blocks, nr_baselines, nr_timesteps);
- idg::Array4D<std::complex<float>> visibilities = idg::get_dummy_visibilities(
- proxy, nr_baselines, nr_timesteps, nr_channels, nr_correlations);
+ aocommon::xt::Span<idg::UVW<float>, 3> uvws =
+ proxy.allocate_span<idg::UVW<float>, 3>(
+ {nr_time_blocks, nr_baselines, nr_timesteps});
+ aocommon::xt::Span<std::complex<float>, 4> visibilities =
+ idg::get_dummy_visibilities(proxy, nr_baselines, nr_timesteps,
+ nr_channels, nr_correlations);
unsigned int bl_offset = 0;
// Vector of plans
@@ -491,8 +496,9 @@ void run_master() {
unsigned int time_offset = t * nr_timesteps;
// Get UVW coordinates for current cycle
- idg::Array2D<idg::UVW<float>> uvw(uvws.data(t, 0, 0), nr_baselines,
- nr_timesteps);
+ const std::array<size_t, 2> uvw_shape{nr_baselines, nr_timesteps};
+ auto uvw = aocommon::xt::CreateSpan(&uvws(t, 0, 0), uvw_shape);
+
if (init) {
runtimes_init[t] -= omp_get_wtime();
data.get_uvw(uvw, bl_offset, time_offset, integration_time);
@@ -638,22 +644,23 @@ void run_worker() {
ProxyType proxy;
// Allocate and initialize static data structures
- idg::Array4D<idg::Matrix2x2<std::complex<float>>> aterms =
+ aocommon::xt::Span<idg::Matrix2x2<std::complex<float>>, 4> aterms =
idg::get_identity_aterms(proxy, nr_timeslots, nr_stations, subgrid_size,
subgrid_size);
- idg::Array1D<unsigned int> aterm_offsets =
+ aocommon::xt::Span<unsigned int, 1> aterm_offsets =
idg::get_example_aterm_offsets(proxy, nr_timeslots, nr_timesteps);
- idg::Array2D<float> spheroidal =
+ aocommon::xt::Span<float, 2> spheroidal =
idg::get_example_spheroidal(proxy, subgrid_size, subgrid_size);
- auto grid =
+ std::shared_ptr<idg::Grid> grid =
proxy.allocate_grid(nr_w_layers, nr_correlations, grid_size, grid_size);
- idg::Array1D<float> shift = idg::get_zero_shift();
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines =
+ std::array<float, 2> shift{0.0f, 0.0f};
+ aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1> baselines =
idg::get_example_baselines(proxy, nr_stations, nr_baselines);
// Receive frequencies
- idg::Array1D<float> frequencies = proxy.allocate_array1d<float>(nr_channels);
- receive_array(0, frequencies);
+ aocommon::xt::Span<float, 1> frequencies =
+ proxy.allocate_span<float, 1>({nr_channels});
+ receive_tensor(0, frequencies);
// Receive data
idg::Data data =
@@ -667,10 +674,12 @@ void run_worker() {
// Buffers for input data
unsigned int nr_time_blocks =
std::ceil((float)total_nr_timesteps / nr_timesteps);
- idg::Array3D<idg::UVW<float>> uvws = proxy.allocate_array3d<idg::UVW<float>>(
- nr_time_blocks, nr_baselines, nr_timesteps);
- idg::Array4D<std::complex<float>> visibilities = idg::get_dummy_visibilities(
- proxy, nr_baselines, nr_timesteps, nr_channels, nr_correlations);
+ aocommon::xt::Span<idg::UVW<float>, 3> uvws =
+ proxy.allocate_span<idg::UVW<float>, 3>(
+ {nr_time_blocks, nr_baselines, nr_timesteps});
+ aocommon::xt::Span<std::complex<float>, 4> visibilities =
+ idg::get_dummy_visibilities(proxy, nr_baselines, nr_timesteps,
+ nr_channels, nr_correlations);
unsigned int bl_offset = 0;
// Vector of plans
@@ -710,8 +719,10 @@ void run_worker() {
unsigned int time_offset = t * nr_timesteps;
// Get UVW coordinates for current cycle
- idg::Array2D<idg::UVW<float>> uvw(uvws.data(t, 0, 0), nr_baselines,
- nr_timesteps);
+ const std::array<size_t, 2> uvw_shape{nr_baselines, nr_timesteps};
+ aocommon::xt::Span<idg::UVW<float>, 2> uvw =
+ aocommon::xt::CreateSpan(&uvws(t, 0, 0), uvw_shape);
+
if (init) {
data.get_uvw(uvw, bl_offset, time_offset, integration_time);
}
diff --git a/idg-bin/examples/cxx/plan/main.cpp b/idg-bin/examples/cxx/plan/main.cpp
index fd788d6d156b81e6140d6f9c602fc338509bfcec..352e560a7163aa617c3e2e652ab836e97c71a8fa 100644
--- a/idg-bin/examples/cxx/plan/main.cpp
+++ b/idg-bin/examples/cxx/plan/main.cpp
@@ -163,22 +163,9 @@ int main(int argc, char** argv) {
std::clog << ">>> Initialize data structures" << std::endl;
// Initialize shift
- idg::Array1D<float> shift = idg::get_zero_shift();
-
- // Initialize frequency data
- idg::Array1D<float> frequencies(nr_channels);
- data.get_frequencies(frequencies, image_size);
-
- // Initalize UVW coordiantes
- idg::Array2D<idg::UVW<float>> uvw(nr_baselines, nr_timesteps);
- data.get_uvw(uvw);
-
- // Initialize metadata
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines =
- idg::get_example_baselines(nr_stations, nr_baselines);
- idg::Array1D<unsigned int> aterm_offsets =
- idg::get_example_aterm_offsets(nr_timeslots, nr_timesteps);
+ std::array<float, 2> shift{0.0f, 0.0f};
+ // Initialize proxy
unsigned int nr_correlations = 4;
idg::proxy::cpu::Optimized proxy;
std::shared_ptr<idg::Grid> grid(
@@ -188,6 +175,21 @@ int main(int argc, char** argv) {
proxy.init_cache(subgrid_size, cell_size, w_step, shift);
+ // Initalize UVW coordinates
+ aocommon::xt::Span<idg::UVW<float>, 2> uvw =
+ proxy.allocate_span<idg::UVW<float>, 2>({nr_baselines, nr_timesteps});
+ data.get_uvw(uvw);
+
+ // Initialize frequency data
+ auto frequencies = proxy.allocate_span<float, 1>({nr_channels});
+ data.get_frequencies(frequencies, image_size);
+
+ // Initialize metadata
+ aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1> baselines =
+ idg::get_example_baselines(proxy, nr_stations, nr_baselines);
+ aocommon::xt::Span<unsigned int, 1> aterm_offsets =
+ idg::get_example_aterm_offsets(proxy, nr_timeslots, nr_timesteps);
+
// Create plan
std::clog << ">>> Create plan" << std::endl;
idg::Plan::Options options;
diff --git a/idg-bin/tests/cxx/CPU/Reference/main.cpp b/idg-bin/tests/cxx/CPU/Reference/main.cpp
index c6bb03134459c1d620ca2d0542f20402a47c45f3..ab9d26efafa0b3f29f66867dd51d7c486f5fdcf9 100644
--- a/idg-bin/tests/cxx/CPU/Reference/main.cpp
+++ b/idg-bin/tests/cxx/CPU/Reference/main.cpp
@@ -1,6 +1,8 @@
// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
+#include <xtensor/xcomplex.hpp>
+
#include "idg-cpu.h" // Reference proxy
#include "idg-util.h" // Data init routines
@@ -50,31 +52,39 @@ int test01() {
// error tolerance, which might need to be adjusted if parameters are changed
float tol = 0.1f;
+ // Initialize proxy
+ clog << ">>> Initialize proxy" << endl;
+ idg::proxy::cpu::Reference proxy;
+
+ // Initialize frequency data
+ aocommon::xt::Span<float, 1> frequencies =
+ proxy.allocate_span<float, 1>({nr_channels});
+ data.get_frequencies(frequencies, image_size);
+
// Allocate and initialize data structures
clog << ">>> Initialize data structures" << endl;
- idg::Array1D<float> frequencies(nr_channels);
- data.get_frequencies(frequencies, image_size);
- idg::Array2D<idg::UVW<float>> uvw = data.get_uvw(nr_baselines, nr_timesteps);
- idg::Array4D<std::complex<float>> visibilities =
- idg::get_example_visibilities(uvw, frequencies, image_size, grid_size,
- nr_correlations);
- idg::Array4D<std::complex<float>> visibilities_ref =
- idg::get_example_visibilities(uvw, frequencies, image_size, grid_size,
- nr_correlations);
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines =
- idg::get_example_baselines(nr_stations, nr_baselines);
+ auto uvw =
+ proxy.allocate_span<idg::UVW<float>, 2>({nr_baselines, nr_timesteps});
+ aocommon::xt::Span<std::complex<float>, 4> visibilities =
+ idg::get_example_visibilities(proxy, uvw, frequencies, image_size,
+ nr_correlations, grid_size);
+ aocommon::xt::Span<std::complex<float>, 4> visibilities_ref =
+ idg::get_example_visibilities(proxy, uvw, frequencies, image_size,
+ nr_correlations, grid_size);
+ aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1> baselines =
+ idg::get_example_baselines(proxy, nr_stations, nr_baselines);
auto grid = std::make_shared<idg::Grid>(nr_w_layers, nr_correlations,
grid_size, grid_size);
auto grid_ref = std::make_shared<idg::Grid>(nr_w_layers, nr_correlations,
grid_size, grid_size);
- idg::Array4D<idg::Matrix2x2<std::complex<float>>> aterms =
- idg::get_identity_aterms(nr_timeslots, nr_stations, subgrid_size,
+ aocommon::xt::Span<idg::Matrix2x2<std::complex<float>>, 4> aterms =
+ idg::get_identity_aterms(proxy, nr_timeslots, nr_stations, subgrid_size,
subgrid_size);
- idg::Array1D<unsigned int> aterm_offsets =
- idg::get_example_aterm_offsets(nr_timeslots, nr_timesteps);
- idg::Array2D<float> spheroidal =
- idg::get_identity_spheroidal(subgrid_size, subgrid_size);
- idg::Array1D<float> shift = idg::get_zero_shift();
+ aocommon::xt::Span<unsigned int, 1> aterm_offsets =
+ idg::get_example_aterm_offsets(proxy, nr_timeslots, nr_timesteps);
+ aocommon::xt::Span<float, 2> spheroidal =
+ idg::get_identity_spheroidal(proxy, subgrid_size, subgrid_size);
+ std::array<float, 2> shift{0.0f, 0.0f};
clog << endl;
// Set w-terms to zero
@@ -96,15 +106,11 @@ int test01() {
(*grid_ref)(0, 1, location_y, location_x) = amplitude;
(*grid_ref)(0, 2, location_y, location_x) = amplitude;
(*grid_ref)(0, 3, location_y, location_x) = amplitude;
- visibilities_ref.zero();
+ visibilities_ref.fill(std::complex<float>(0.0f, 0.0f));
add_pt_src(visibilities_ref, uvw, frequencies, image_size, grid_size,
offset_x, offset_y, amplitude);
clog << endl;
- // Initialize proxy
- clog << ">>> Initialize proxy" << endl;
- idg::proxy::cpu::Reference proxy;
-
// Set grid
proxy.set_grid(grid);
float w_step = 0.0;
diff --git a/idg-bin/tests/cxx/CUDA/wtiling/backward.cpp b/idg-bin/tests/cxx/CUDA/wtiling/backward.cpp
index 25eec893c8476f94d6c04cced289651d5be92337..ad8957444d40d61c6abf2e95c8e5b7aa9aedf866 100644
--- a/idg-bin/tests/cxx/CUDA/wtiling/backward.cpp
+++ b/idg-bin/tests/cxx/CUDA/wtiling/backward.cpp
@@ -141,7 +141,8 @@ void subgrids_from_wtiles(const long nr_subgrids, const int nr_polarizations,
void wtiles_from_grid(int nr_tiles, int nr_polarizations, int grid_size,
int tile_size, int padded_tile_size, int* tile_ids,
idg::Coordinate* tile_coordinates,
- std::complex<float>* tiles, std::complex<float>* grid) {
+ xt::xtensor<std::complex<float>, 4>& tiles,
+ std::complex<float>* grid) {
// Extract tiles from grid
for (int i = 0; i < nr_tiles; i++) {
idg::Coordinate& coordinate = tile_coordinates[i];
@@ -161,12 +162,8 @@ void wtiles_from_grid(int nr_tiles, int nr_polarizations, int grid_size,
for (int pol = 0; pol < nr_polarizations; pol++) {
const int index_pol_transposed[4] = {0, 2, 1, 3};
unsigned int pol_src = index_pol_transposed[pol];
- unsigned int pol_dst = pol;
unsigned long src_idx = index_grid_3d(grid_size, pol_src, y, x);
- unsigned long dst_idx =
- index_grid_4d(nr_polarizations, padded_tile_size, i, pol_dst,
- (y - y0), (x - x0));
- tiles[dst_idx] = grid[src_idx];
+ tiles(i, pol, y - y0, x - x0) = grid[src_idx];
}
}
}
@@ -200,17 +197,22 @@ int main(int argc, char* argv[]) {
float cell_size = image_size / grid_size;
// Initialize data
- idg::Array1D<float> frequencies(nr_channels);
- data.get_frequencies(frequencies, image_size);
- idg::Array2D<idg::UVW<float>> uvw(nr_baselines, nr_timesteps);
- data.get_uvw(uvw);
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines =
+ const std::array<size_t, 1> frequencies_shape{nr_channels};
+ xt::xtensor<float, 1> frequencies(frequencies_shape);
+ aocommon::xt::Span<float, 1> frequencies_span =
+ aocommon::xt::CreateSpan(frequencies);
+ data.get_frequencies(frequencies_span, image_size);
+ const std::array<size_t, 2> uvw_shape{nr_baselines, nr_timesteps};
+ xt::xtensor<idg::UVW<float>, 2> uvw(uvw_shape);
+ auto uvw_span = aocommon::xt::CreateSpan(uvw);
+ data.get_uvw(uvw_span);
+ xt::xtensor<std::pair<unsigned int, unsigned int>, 1> baselines =
idg::get_example_baselines(nr_stations, nr_baselines);
- idg::Array1D<unsigned int> aterm_offsets =
- idg::get_example_aterm_offsets(nr_timeslots, nr_timesteps);
- idg::Array1D<float> shift = idg::get_zero_shift();
- shift(0) = 10.1f;
- shift(0) = 20.2f;
+ auto baselines_span = aocommon::xt::CreateSpan(baselines);
+ xt::xtensor<unsigned int, 1> aterm_offsets({nr_timeslots + 1}, 0);
+ auto aterm_offsets_span = aocommon::xt::CreateSpan(aterm_offsets);
+ idg::init_example_aterm_offsets(aterm_offsets_span, nr_timesteps);
+ std::array<float, 2> shift{10.1f, 20.2f};
// Set w-terms to zero
for (unsigned bl = 0; bl < nr_baselines; bl++) {
@@ -224,7 +226,8 @@ int main(int argc, char* argv[]) {
idg::Plan::Options options;
options.plan_strict = true;
idg::Plan plan(kernel_size, subgrid_size, grid_size, cell_size, shift,
- frequencies, uvw, baselines, aterm_offsets, wtiles, options);
+ frequencies_span, uvw_span, baselines_span, aterm_offsets_span,
+ wtiles, options);
int nr_subgrids = plan.get_nr_subgrids();
// Get W-Tiling paramters
@@ -305,7 +308,7 @@ int main(int argc, char* argv[]) {
h_tiles.size());
// Init shift
- stream.memcpyHtoDAsync(d_shift, shift.data(), shift.bytes());
+ stream.memcpyHtoDAsync(d_shift, shift.data(), sizeof_shift);
// Start tests
std::cout << std::endl;
@@ -369,13 +372,13 @@ int main(int argc, char* argv[]) {
stream.synchronize();
// Run splitter_wtiles_from_grid on host
- idg::Array4D<std::complex<float>> padded_tiles(
- max_nr_tiles, nr_polarizations, padded_tile_size, padded_tile_size);
- padded_tiles.zero();
+ xt::xtensor<std::complex<float>, 4> padded_tiles(
+ {max_nr_tiles, nr_polarizations, padded_tile_size, padded_tile_size},
+ std::complex<float>(0.0f, 0.0f));
wtiles_from_grid(nr_tiles, nr_polarizations, grid_size,
tile_size - subgrid_size, padded_tile_size,
static_cast<int*>(h_padded_tile_ids.data()),
- tile_coordinates.data(), padded_tiles.data(),
+ tile_coordinates.data(), padded_tiles,
static_cast<std::complex<float>*>(u_grid.data()));
n = nr_tiles * nr_polarizations * padded_tile_size * padded_tile_size;
diff --git a/idg-bin/tests/cxx/CUDA/wtiling/forward.cpp b/idg-bin/tests/cxx/CUDA/wtiling/forward.cpp
index 9433eef08f5e3975ec31f693ae75cc94e54e9c8a..095203d509436f2a3191eafec9bfacfc4b93f709 100644
--- a/idg-bin/tests/cxx/CUDA/wtiling/forward.cpp
+++ b/idg-bin/tests/cxx/CUDA/wtiling/forward.cpp
@@ -217,7 +217,8 @@ void subgrids_to_wtiles(const long nr_subgrids, const int nr_polarizations,
void wtiles_to_grid(int nr_tiles, int nr_polarizations, int grid_size,
int tile_size, int padded_tile_size, int* tile_ids,
idg::Coordinate* tile_coordinates,
- std::complex<float>* tiles, std::complex<float>* grid) {
+ std::complex<float>* tiles,
+ xt::xtensor<std::complex<float>, 3>& grid) {
for (int i = 0; i < nr_tiles; i++) {
idg::Coordinate& coordinate = tile_coordinates[i];
@@ -236,11 +237,10 @@ void wtiles_to_grid(int nr_tiles, int nr_polarizations, int grid_size,
const int index_pol_transposed[nr_polarizations] = {0, 2, 1, 3};
unsigned int pol_src = pol;
unsigned int pol_dst = index_pol_transposed[pol];
- unsigned long dst_idx = index_grid_3d(grid_size, pol_src, y, x);
unsigned long src_idx =
index_grid_4d(nr_polarizations, padded_tile_size, i, pol_dst,
(y - y0), (x - x0));
- grid[dst_idx] += tiles[src_idx];
+ grid(pol_src, y, x) += tiles[src_idx];
} // end for pol
} // end for x
} // end for y
@@ -274,17 +274,21 @@ int main(int argc, char* argv[]) {
float w_step = 4.0 / (image_size * image_size);
// Initialize data
- idg::Array1D<float> frequencies(nr_channels);
- data.get_frequencies(frequencies, image_size);
- idg::Array2D<idg::UVW<float>> uvw(nr_baselines, nr_timesteps);
- data.get_uvw(uvw);
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines =
+ const std::array<size_t, 1> frequencies_shape{nr_channels};
+ xt::xtensor<float, 1> frequencies(frequencies_shape);
+ aocommon::xt::Span<float, 1> frequencies_span =
+ aocommon::xt::CreateSpan(frequencies);
+ data.get_frequencies(frequencies_span, image_size);
+ const std::array<size_t, 2> uvw_shape{nr_baselines, nr_timesteps};
+ xt::xtensor<idg::UVW<float>, 2> uvw(uvw_shape);
+ auto uvw_span = aocommon::xt::CreateSpan(uvw);
+ data.get_uvw(uvw_span);
+ xt::xtensor<std::pair<unsigned int, unsigned int>, 1> baselines =
idg::get_example_baselines(nr_stations, nr_baselines);
- idg::Array1D<unsigned int> aterm_offsets =
- idg::get_example_aterm_offsets(nr_timeslots, nr_timesteps);
- idg::Array1D<float> shift = idg::get_zero_shift();
- shift(0) = 10.1f;
- shift(0) = 20.2f;
+ auto baselines_span = aocommon::xt::CreateSpan(baselines);
+ xt::xtensor<unsigned int, 1> aterm_offsets({nr_timeslots + 1}, 0);
+ auto aterm_offsets_span = aocommon::xt::CreateSpan(aterm_offsets);
+ std::array<float, 2> shift{10.1f, 20.2f};
// Set w-terms to zero
for (unsigned bl = 0; bl < nr_baselines; bl++) {
@@ -298,7 +302,8 @@ int main(int argc, char* argv[]) {
idg::Plan::Options options;
options.plan_strict = true;
idg::Plan plan(kernel_size, subgrid_size, grid_size, cell_size, shift,
- frequencies, uvw, baselines, aterm_offsets, wtiles, options);
+ frequencies_span, uvw_span, baselines_span, aterm_offsets_span,
+ wtiles, options);
int nr_subgrids = plan.get_nr_subgrids();
// Get W-Tiling paramters
@@ -379,7 +384,7 @@ int main(int argc, char* argv[]) {
h_tiles.size());
// Init shift
- stream.memcpyHtoDAsync(d_shift, shift.data(), shift.bytes());
+ stream.memcpyHtoDAsync(d_shift, shift.data(), sizeof_shift);
// Start tests
std::cout << std::endl;
@@ -549,14 +554,14 @@ int main(int argc, char* argv[]) {
stream.synchronize();
// Run adder_wtiles_to_grid on host
- idg::Array3D<std::complex<float>> grid(nr_polarizations, grid_size,
- grid_size);
- grid.zero();
+ xt::xtensor<std::complex<float>, 3> grid(
+ {nr_polarizations, grid_size, grid_size},
+ std::complex<float>(0.0f, 0.0f));
wtiles_to_grid(
nr_tiles, nr_polarizations, grid_size, tile_size - subgrid_size,
padded_tile_size, static_cast<int*>(h_padded_tile_ids.data()),
tile_coordinates.data(),
- static_cast<std::complex<float>*>(h_padded_tiles.data()), grid.data());
+ static_cast<std::complex<float>*>(h_padded_tiles.data()), grid);
n = nr_polarizations * grid_size * grid_size;
accuracy = compare_arrays(n, static_cast<std::complex<float>*>(u_grid.data()),
@@ -632,12 +637,12 @@ int main(int argc, char* argv[]) {
}
// Create reference grid
- grid.zero();
+ grid.fill(std::complex<float>(0.0f, 0.0f));
wtiles_to_grid(
nr_tiles, nr_polarizations, grid_size, tile_size - subgrid_size,
padded_tile_size, static_cast<int*>(h_padded_tile_ids.data()),
tile_coordinates.data(),
- static_cast<std::complex<float>*>(h_padded_tiles.data()), grid.data());
+ static_cast<std::complex<float>*>(h_padded_tiles.data()), grid);
n = nr_polarizations * grid_size * grid_size;
accuracy = compare_arrays(n, static_cast<std::complex<float>*>(u_grid.data()),
diff --git a/idg-bin/tests/cxx/common/common.cpp b/idg-bin/tests/cxx/common/common.cpp
index 87d301098d6e5ffeb9e0757addcf7b2b1073b695..28f9187c570e87412779e212507fc0791fffbed1 100644
--- a/idg-bin/tests/cxx/common/common.cpp
+++ b/idg-bin/tests/cxx/common/common.cpp
@@ -145,30 +145,30 @@ int compare(idg::proxy::Proxy& proxy1, idg::proxy::Proxy& proxy2, float tol) {
image_size, grid_size, subgrid_size, kernel_size);
std::clog << ">>> Initialize data structures" << std::endl;
- idg::Array1D<float> frequencies = proxy2.allocate_array1d<float>(nr_channels);
-
+ aocommon::xt::Span<float, 1> frequencies =
+ proxy2.allocate_span<float, 1>({nr_channels});
data.get_frequencies(frequencies, image_size);
- idg::Array2D<idg::UVW<float>> uvw =
- proxy2.allocate_array2d<idg::UVW<float>>(nr_baselines, nr_timesteps);
+ aocommon::xt::Span<idg::UVW<float>, 2> uvw =
+ proxy2.allocate_span<idg::UVW<float>, 2>({nr_baselines, nr_timesteps});
data.get_uvw(uvw);
- idg::Array4D<std::complex<float>> visibilities = idg::get_dummy_visibilities(
- proxy2, nr_baselines, nr_timesteps, nr_channels, nr_correlations);
- idg::Array4D<std::complex<float>> visibilities_ref =
+ aocommon::xt::Span<std::complex<float>, 4> visibilities =
+ idg::get_dummy_visibilities(proxy2, nr_baselines, nr_timesteps,
+ nr_channels, nr_correlations);
+ aocommon::xt::Span<std::complex<float>, 4> visibilities_ref =
idg::get_dummy_visibilities(proxy1, nr_baselines, nr_timesteps,
nr_channels, nr_correlations);
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines =
+ aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1> baselines =
idg::get_example_baselines(proxy2, nr_stations, nr_baselines);
- idg::Array4D<idg::Matrix2x2<std::complex<float>>> aterms =
+ aocommon::xt::Span<idg::Matrix2x2<std::complex<float>>, 4> aterms =
idg::get_example_aterms(proxy2, nr_timeslots, nr_stations, subgrid_size,
subgrid_size);
- idg::Array1D<unsigned int> aterm_offsets =
+ aocommon::xt::Span<unsigned int, 1> aterm_offsets =
idg::get_example_aterm_offsets(proxy2, nr_timeslots, nr_timesteps);
- idg::Array2D<float> spheroidal =
+ aocommon::xt::Span<float, 2> spheroidal =
idg::get_example_spheroidal(proxy2, subgrid_size, subgrid_size);
- idg::Array1D<float> shift = idg::get_zero_shift();
// Set shift to some random non-zero value
- shift(0) = 0.02; // shift is in radians, 0.02rad is about 1deg
- shift(1) = -0.03;
+ // shift is in radians, 0.02rad is about 1deg
+ std::array<float, 2> shift{0.02f, -0.03f};
auto grid = proxy2.allocate_grid(1, nr_polarizations, grid_size, grid_size);
auto grid_ref =
proxy1.allocate_grid(1, nr_polarizations, grid_size, grid_size);
@@ -237,12 +237,12 @@ int compare(idg::proxy::Proxy& proxy1, idg::proxy::Proxy& proxy2, float tol) {
#if TEST_DEGRIDDING
// Run degridder
std::clog << ">>> Run degridding" << std::endl;
- visibilities.zero();
+ visibilities.fill(std::complex<float>(0.0f, 0.0f));
proxy2.degridding(*plan2, frequencies, visibilities, uvw, baselines, aterms,
aterm_offsets, spheroidal);
std::clog << ">>> Run reference degridding" << std::endl;
- visibilities_ref.zero();
+ visibilities_ref.fill(std::complex<float>(0.0f, 0.0f));
proxy1.degridding(*plan1, frequencies, visibilities_ref, uvw, baselines,
aterms, aterm_offsets, spheroidal);
diff --git a/idg-lib/src/CPU/common/CPU.cpp b/idg-lib/src/CPU/common/CPU.cpp
index b9f9c1afdf3e6e75a537e06667a3326f288c5f08..ef527a248f1cc9bfd32c3cf2908cfa017b87b260 100644
--- a/idg-lib/src/CPU/common/CPU.cpp
+++ b/idg-lib/src/CPU/common/CPU.cpp
@@ -43,10 +43,12 @@ std::unique_ptr<auxiliary::Memory> CPU::allocate_memory(size_t bytes) {
}
std::unique_ptr<Plan> CPU::make_plan(
- const int kernel_size, const Array1D<float>& frequencies,
- const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets, Plan::Options options) {
+ const int kernel_size, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
+ Plan::Options options) {
if (supports_wtiling() && m_cache_state.w_step != 0.0 &&
m_wtiles.get_wtile_buffer_size()) {
options.w_step = m_cache_state.w_step;
@@ -62,7 +64,7 @@ std::unique_ptr<Plan> CPU::make_plan(
}
void CPU::init_cache(int subgrid_size, float cell_size, float w_step,
- const Array1D<float>& shift) {
+ const std::array<float, 2>& shift) {
Proxy::init_cache(subgrid_size, cell_size, w_step, shift);
const int nr_polarizations = m_grid->get_z_dim();
const size_t grid_size = m_grid->get_x_dim();
diff --git a/idg-lib/src/CPU/common/CPU.h b/idg-lib/src/CPU/common/CPU.h
index ccf293575cda46a1888d9aaa7c56f881b27b6da4..1207e51a2dec5ca9f1dd676aa0cc613c15ba91e9 100644
--- a/idg-lib/src/CPU/common/CPU.h
+++ b/idg-lib/src/CPU/common/CPU.h
@@ -37,14 +37,15 @@ class CPU : public Proxy {
std::shared_ptr<kernel::cpu::InstanceCPU> get_kernels() { return m_kernels; }
std::unique_ptr<Plan> make_plan(
- const int kernel_size, const Array1D<float>& frequencies,
- const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets,
+ const int kernel_size, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
Plan::Options options) override;
void init_cache(int subgrid_size, float cell_size, float w_step,
- const Array1D<float>& shift) override;
+ const std::array<float, 2>& shift) override;
std::shared_ptr<Grid> get_final_grid() override;
diff --git a/idg-lib/src/CUDA/Generic/Generic.cpp b/idg-lib/src/CUDA/Generic/Generic.cpp
index 21d292ad6d2bbb16368f18080789bb1630b16661..7bd45cd1e81c10b011ec629669dc4a39328695e9 100644
--- a/idg-lib/src/CUDA/Generic/Generic.cpp
+++ b/idg-lib/src/CUDA/Generic/Generic.cpp
@@ -123,10 +123,12 @@ std::shared_ptr<Grid> Generic::get_final_grid() {
}
std::unique_ptr<Plan> Generic::make_plan(
- const int kernel_size, const Array1D<float>& frequencies,
- const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets, Plan::Options options) {
+ const int kernel_size, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
+ Plan::Options options) {
if (do_supports_wtiling() && !m_disable_wtiling) {
options.w_step = m_cache_state.w_step;
options.nr_w_layers = std::numeric_limits<int>::max();
@@ -141,7 +143,7 @@ std::unique_ptr<Plan> Generic::make_plan(
}
void Generic::init_cache(int subgrid_size, float cell_size, float w_step,
- const Array1D<float>& shift) {
+ const std::array<float, 2>& shift) {
// Initialize cache
Proxy::init_cache(subgrid_size, cell_size, w_step, shift);
diff --git a/idg-lib/src/CUDA/Generic/Generic.h b/idg-lib/src/CUDA/Generic/Generic.h
index e1ece0bc5e62bb3237e80abf8662a3959aaed08d..13db417830537138d45cfd3acd3268090e1dd29d 100644
--- a/idg-lib/src/CUDA/Generic/Generic.h
+++ b/idg-lib/src/CUDA/Generic/Generic.h
@@ -64,14 +64,15 @@ class Generic : public CUDA {
std::shared_ptr<Grid> get_final_grid() override;
virtual std::unique_ptr<Plan> make_plan(
- const int kernel_size, const Array1D<float>& frequencies,
- const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets,
+ const int kernel_size, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
Plan::Options options) override;
void init_cache(int subgrid_size, float cell_size, float w_step,
- const Array1D<float>& shift) override;
+ const std::array<float, 2>& shift) override;
private:
std::unique_ptr<cu::DeviceMemory> d_grid_;
diff --git a/idg-lib/src/CUDA/Generic/routines/Imaging.cpp b/idg-lib/src/CUDA/Generic/routines/Imaging.cpp
index 3490c43b4d4d82294f035f43258f0dc129b50ddb..d5d0ef55f28a6a50d86a2e1c610c7acca1574768 100644
--- a/idg-lib/src/CUDA/Generic/routines/Imaging.cpp
+++ b/idg-lib/src/CUDA/Generic/routines/Imaging.cpp
@@ -242,8 +242,8 @@ void Generic::run_imaging(
// Launch gridder kernel
device.launch_gridder(
time_offset_current, nr_subgrids_current, nr_polarizations, grid_size,
- subgrid_size, image_size, w_step, nr_channels, nr_stations, shift(0),
- shift(1), d_uvw, d_wavenumbers, d_visibilities, d_spheroidal,
+ subgrid_size, image_size, w_step, nr_channels, nr_stations, shift[0],
+ shift[1], d_uvw, d_wavenumbers, d_visibilities, d_spheroidal,
d_aterms, d_aterm_indices, d_metadata, *d_avg_aterm, d_subgrids);
// Launch FFT
@@ -294,8 +294,8 @@ void Generic::run_imaging(
// Launch degridder kernel
device.launch_degridder(
time_offset_current, nr_subgrids_current, nr_polarizations, grid_size,
- subgrid_size, image_size, w_step, nr_channels, nr_stations, shift(0),
- shift(1), d_uvw, d_wavenumbers, d_visibilities, d_spheroidal,
+ subgrid_size, image_size, w_step, nr_channels, nr_stations, shift[0],
+ shift[1], d_uvw, d_wavenumbers, d_visibilities, d_spheroidal,
d_aterms, d_aterm_indices, d_metadata, d_subgrids);
executestream.record(gpuFinished[job_id]);
diff --git a/idg-lib/src/CUDA/common/CUDA.cpp b/idg-lib/src/CUDA/common/CUDA.cpp
index 0e15faee9af8223a2370388cd4d875f743736891..45cf269e7b249e38d8a0b38054b6ffd2ed823184 100644
--- a/idg-lib/src/CUDA/common/CUDA.cpp
+++ b/idg-lib/src/CUDA/common/CUDA.cpp
@@ -32,6 +32,7 @@ CUDA::~CUDA() {
// contexts are free'ed, hence the explicit calls here.
free_unified_grid();
free_buffers_wtiling();
+ free_memory();
}
void CUDA::init_devices() {
diff --git a/idg-lib/src/CUDA/common/CUDA.h b/idg-lib/src/CUDA/common/CUDA.h
index 8a4b9f467ea0e9f0ac7edc3a21db6d69a79bd03d..18d9761583e2e2e655bb05fd4afff2cd5748e1dd 100644
--- a/idg-lib/src/CUDA/common/CUDA.h
+++ b/idg-lib/src/CUDA/common/CUDA.h
@@ -119,31 +119,31 @@ class CUDA : public Proxy {
const int nr_polarizations,
const int subgrid_size, const float image_size,
const float w_step,
- const idg::Array1D<float>& shift,
+ const std::array<float, 2>& shift,
const size_t bytes_free) const;
void run_wtiles_to_grid(unsigned int subgrid_size, float image_size,
- float w_step, const Array1D<float>& shift,
+ float w_step, const std::array<float, 2>& shift,
WTileUpdateInfo& wtile_flush_info);
void run_subgrids_to_wtiles(unsigned int nr_polarizations,
unsigned int subgrid_offset,
unsigned int nr_subgrids,
unsigned int subgrid_size, float image_size,
- float w_step, const Array1D<float>& shift,
+ float w_step, const std::array<float, 2>& shift,
WTileUpdateSet& wtile_flush_set,
cu::DeviceMemory& d_subgrids,
cu::DeviceMemory& d_metadata);
void run_wtiles_from_grid(unsigned int subgrid_size, float image_size,
- float w_step, const Array1D<float>& shift,
+ float w_step, const std::array<float, 2>& shift,
WTileUpdateInfo& wtile_initialize_info);
void run_subgrids_from_wtiles(unsigned int nr_polarizations,
unsigned int subgrid_offset,
unsigned int nr_subgrids,
unsigned int subgrid_size, float image_size,
- float w_step, const Array1D<float>& shift,
+ float w_step, const std::array<float, 2>& shift,
WTileUpdateSet& wtile_initialize_set,
cu::DeviceMemory& d_subgrids,
cu::DeviceMemory& d_metadata);
diff --git a/idg-lib/src/CUDA/common/routines/WTiling.cpp b/idg-lib/src/CUDA/common/routines/WTiling.cpp
index f5fd663352b7f1e20d0f543f5de4a0130cd30993..9454ecff2e8882f69ae1b35b60b89948a9bdbcb3 100644
--- a/idg-lib/src/CUDA/common/routines/WTiling.cpp
+++ b/idg-lib/src/CUDA/common/routines/WTiling.cpp
@@ -55,16 +55,15 @@ size_t CUDA::bytes_required_wtiling(const WTileUpdateSet& wtile_set,
const int nr_polarizations,
const int subgrid_size,
const float image_size, const float w_step,
- const idg::Array1D<float>& shift,
+ const std::array<float, 2>& shift,
const size_t bytes_free) const {
if (wtile_set.empty()) {
return 0;
}
// Compute the maximum padded w-tile size for any w-tile in the set
- int w_padded_tile_size =
- compute_w_padded_tile_size_max(wtile_set, m_tile_size, subgrid_size,
- image_size, w_step, shift(0), shift(1));
+ int w_padded_tile_size = compute_w_padded_tile_size_max(
+ wtile_set, m_tile_size, subgrid_size, image_size, w_step, shift);
// Compute the memory required for such a padded w-tile
size_t sizeof_w_padded_tile = w_padded_tile_size * w_padded_tile_size *
@@ -159,7 +158,7 @@ void CUDA::init_buffers_wtiling(unsigned int subgrid_size) {
}
void CUDA::run_wtiles_to_grid(unsigned int subgrid_size, float image_size,
- float w_step, const Array1D<float>& shift,
+ float w_step, const std::array<float, 2>& shift,
WTileUpdateInfo& wtile_flush_info) {
// Load grid parameters
unsigned int nr_polarizations = m_grid->get_z_dim();
@@ -192,7 +191,7 @@ void CUDA::run_wtiles_to_grid(unsigned int subgrid_size, float image_size,
// Compute w_padded_tile_size for all tiles
const float image_size_shift =
- image_size + 2 * std::max(std::abs(shift(0)), std::abs(shift(1)));
+ image_size + 2 * std::max(std::abs(shift[0]), std::abs(shift[1]));
std::vector<int> w_padded_tile_sizes = compute_w_padded_tile_sizes(
tile_coordinates.data(), nr_tiles, w_step, image_size, image_size_shift,
padded_tile_size);
@@ -225,8 +224,9 @@ void CUDA::run_wtiles_to_grid(unsigned int subgrid_size, float image_size,
}
// Copy shift to device
- cu::DeviceMemory d_shift(context, shift.bytes());
- executestream.memcpyHtoDAsync(d_shift, shift.data(), shift.bytes());
+ const size_t sizeof_shift = shift.size() * sizeof(float);
+ cu::DeviceMemory d_shift(context, sizeof_shift);
+ executestream.memcpyHtoDAsync(d_shift, shift.data(), sizeof_shift);
// FFT plan
std::unique_ptr<cufft::C2C_2D> fft;
@@ -386,7 +386,7 @@ void CUDA::run_wtiles_to_grid(unsigned int subgrid_size, float image_size,
void CUDA::run_subgrids_to_wtiles(
unsigned nr_polarizations, unsigned int subgrid_offset,
unsigned int nr_subgrids, unsigned int subgrid_size, float image_size,
- float w_step, const idg::Array1D<float>& shift,
+ float w_step, const std::array<float, 2>& shift,
WTileUpdateSet& wtile_flush_set, cu::DeviceMemory& d_subgrids,
cu::DeviceMemory& d_metadata) {
// Load CUDA objects
@@ -447,7 +447,7 @@ void CUDA::run_subgrids_to_wtiles(
}
void CUDA::run_wtiles_from_grid(unsigned int subgrid_size, float image_size,
- float w_step, const Array1D<float>& shift,
+ float w_step, const std::array<float, 2>& shift,
WTileUpdateInfo& wtile_initialize_info) {
// Load grid parameters
unsigned int nr_polarizations = m_grid->get_z_dim();
@@ -480,7 +480,7 @@ void CUDA::run_wtiles_from_grid(unsigned int subgrid_size, float image_size,
// Compute w_padded_tile_size for all tiles
const float image_size_shift =
- image_size + 2 * std::max(std::abs(shift(0)), std::abs(shift(1)));
+ image_size + 2 * std::max(std::abs(shift[0]), std::abs(shift[1]));
std::vector<int> w_padded_tile_sizes = compute_w_padded_tile_sizes(
tile_coordinates.data(), nr_tiles, w_step, image_size, image_size_shift,
padded_tile_size);
@@ -517,8 +517,9 @@ void CUDA::run_wtiles_from_grid(unsigned int subgrid_size, float image_size,
sizeof_tile_ids);
// Copy shift to device
- cu::DeviceMemory d_shift(context, shift.bytes());
- executestream.memcpyHtoDAsync(d_shift, shift.data(), shift.bytes());
+ const size_t sizeof_shift = shift.size() * sizeof(float);
+ cu::DeviceMemory d_shift(context, sizeof_shift);
+ executestream.memcpyHtoDAsync(d_shift, shift.data(), sizeof_shift);
// FFT plan
std::unique_ptr<cufft::C2C_2D> fft;
@@ -675,14 +676,12 @@ void CUDA::run_wtiles_from_grid(unsigned int subgrid_size, float image_size,
} // end for tile_offset
}
-void CUDA::run_subgrids_from_wtiles(unsigned int nr_polarizations,
- unsigned int subgrid_offset,
- unsigned int nr_subgrids,
- unsigned int subgrid_size, float image_size,
- float w_step, const Array1D<float>& shift,
- WTileUpdateSet& wtile_initialize_set,
- cu::DeviceMemory& d_subgrids,
- cu::DeviceMemory& d_metadata) {
+void CUDA::run_subgrids_from_wtiles(
+ unsigned int nr_polarizations, unsigned int subgrid_offset,
+ unsigned int nr_subgrids, unsigned int subgrid_size, float image_size,
+ float w_step, const std::array<float, 2>& shift,
+ WTileUpdateSet& wtile_initialize_set, cu::DeviceMemory& d_subgrids,
+ cu::DeviceMemory& d_metadata) {
// Load CUDA objects
kernel::cuda::InstanceCUDA& device = get_device(0);
cu::Stream& stream = device.get_execute_stream();
@@ -749,7 +748,7 @@ void CUDA::flush_wtiles() {
float image_size = grid_size * cell_size;
int subgrid_size = m_cache_state.subgrid_size;
float w_step = m_cache_state.w_step;
- const Array1D<float>& shift = m_cache_state.shift;
+ const std::array<float, 2>& shift = m_cache_state.shift;
// Get all the remaining wtiles
WTileUpdateInfo wtile_flush_info = m_wtiles.clear();
diff --git a/idg-lib/src/Hybrid/CUDA/GenericOptimized/GenericOptimized.cpp b/idg-lib/src/Hybrid/CUDA/GenericOptimized/GenericOptimized.cpp
index cd236fbf2020a2c762c00561c2143e1d14462271..5f60767ab66a7a2275a795d03fbd843c35880203 100644
--- a/idg-lib/src/Hybrid/CUDA/GenericOptimized/GenericOptimized.cpp
+++ b/idg-lib/src/Hybrid/CUDA/GenericOptimized/GenericOptimized.cpp
@@ -63,10 +63,12 @@ GenericOptimized::~GenericOptimized() {
* Plan
*/
std::unique_ptr<Plan> GenericOptimized::make_plan(
- const int kernel_size, const Array1D<float>& frequencies,
- const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets, Plan::Options options) {
+ const int kernel_size, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
+ Plan::Options options) {
if (!m_disable_wtiling && !m_disable_wtiling_gpu) {
options.w_step = m_cache_state.w_step;
options.nr_w_layers = INT_MAX;
@@ -165,7 +167,8 @@ std::shared_ptr<Grid> GenericOptimized::get_final_grid() {
* Cache
*/
void GenericOptimized::init_cache(int subgrid_size, float cell_size,
- float w_step, const Array1D<float>& shift) {
+ float w_step,
+ const std::array<float, 2>& shift) {
// Initialize cache
Proxy::init_cache(subgrid_size, cell_size, w_step, shift);
diff --git a/idg-lib/src/Hybrid/CUDA/GenericOptimized/GenericOptimized.h b/idg-lib/src/Hybrid/CUDA/GenericOptimized/GenericOptimized.h
index 404e57d932f6b3dcb2cb6284c9074ce02e0f78f8..0d4e9ac776a12037ec2bd13faa1884c598f228fe 100644
--- a/idg-lib/src/Hybrid/CUDA/GenericOptimized/GenericOptimized.h
+++ b/idg-lib/src/Hybrid/CUDA/GenericOptimized/GenericOptimized.h
@@ -44,13 +44,14 @@ class GenericOptimized : public cuda::CUDA {
std::shared_ptr<Grid> get_final_grid() override;
void init_cache(int subgrid_size, float cell_size, float w_step,
- const Array1D<float>& shift) override;
+ const std::array<float, 2>& shift) override;
std::unique_ptr<Plan> make_plan(
- const int kernel_size, const Array1D<float>& frequencies,
- const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets,
+ const int kernel_size, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
Plan::Options options) override;
private:
diff --git a/idg-lib/src/Hybrid/CUDA/GenericOptimized/routines/Imaging.cpp b/idg-lib/src/Hybrid/CUDA/GenericOptimized/routines/Imaging.cpp
index 219e6586eb6d74b72fdec54ce77369318f881f8f..63e43aec8cb3177727d66d8597e8a183b220907a 100644
--- a/idg-lib/src/Hybrid/CUDA/GenericOptimized/routines/Imaging.cpp
+++ b/idg-lib/src/Hybrid/CUDA/GenericOptimized/routines/Imaging.cpp
@@ -252,8 +252,8 @@ void GenericOptimized::run_imaging(
// Launch gridder kernel
device.launch_gridder(
time_offset_current, nr_subgrids_current, nr_polarizations, grid_size,
- subgrid_size, image_size, w_step, nr_channels, nr_stations, shift(0),
- shift(1), d_uvw, d_wavenumbers, d_visibilities, d_spheroidal,
+ subgrid_size, image_size, w_step, nr_channels, nr_stations, shift[0],
+ shift[1], d_uvw, d_wavenumbers, d_visibilities, d_spheroidal,
d_aterms, d_aterm_indices, d_metadata, *d_avg_aterm, d_subgrids);
// Launch FFT
@@ -352,8 +352,8 @@ void GenericOptimized::run_imaging(
// Launch degridder kernel
device.launch_degridder(
time_offset_current, nr_subgrids_current, nr_polarizations, grid_size,
- subgrid_size, image_size, w_step, nr_channels, nr_stations, shift(0),
- shift(1), d_uvw, d_wavenumbers, d_visibilities, d_spheroidal,
+ subgrid_size, image_size, w_step, nr_channels, nr_stations, shift[0],
+ shift[1], d_uvw, d_wavenumbers, d_visibilities, d_spheroidal,
d_aterms, d_aterm_indices, d_metadata, d_subgrids);
executestream.record(gpuFinished[job_id]);
diff --git a/idg-lib/src/common/ArrayTypes.h b/idg-lib/src/common/ArrayTypes.h
index 283cb80bcbf579aad24f6cb3bc2432cd4c52afa2..de64005446ef28ba7c12932b138337a26e77ff59 100644
--- a/idg-lib/src/common/ArrayTypes.h
+++ b/idg-lib/src/common/ArrayTypes.h
@@ -8,6 +8,7 @@
#include <cassert>
#include <omp.h>
+#include <aocommon/xt/span.h>
#include "auxiliary.h"
#include "Types.h"
@@ -47,6 +48,20 @@ class Array1D {
Array1D(T* data, size_t size)
: m_x_dim(size), m_memory(nullptr), m_buffer(data) {}
+ /**
+ * @brief Construct a new Array1D object from a Span
+ *
+ * @param data The Span that the Array will use.
+ * Must be avaible during the life time of the Array
+ * The user is responsible for allocation and deallocation.
+ */
+ Array1D(const aocommon::xt::Span<T, 1>& data)
+ : m_x_dim(data.size()),
+ m_memory(nullptr),
+ // This const_cast simplifies the transition to XTensor. It will
+ // be removed at the end of the transition, together with Array1D.
+ m_buffer(const_cast<T*>(data.data())) {}
+
/**
* @brief Construct a new Array 1 D object from a Memory object
*
@@ -263,6 +278,21 @@ class Array2D : public Array1D<T> {
Array2D(T* data, size_t y_dim, size_t x_dim)
: Array1D<T>(data, y_dim * x_dim), m_y_dim(y_dim), m_x_dim(x_dim) {}
+ /**
+ * @brief Construct a new Array2D object from a Span
+ *
+ * @param data The Span that the Array will use.
+ * Must be avaible during the life time of the Array
+ * The user is responsible for allocation and deallocation.
+ */
+ Array2D(const aocommon::xt::Span<T, 2>& data)
+ : m_y_dim(data.shape(0)), m_x_dim(data.shape(1)) {
+ this->m_memory = nullptr;
+ // This const_cast simplifies the transition to XTensor. It will
+ // be removed at the end of the transition, together with Array2D.
+ this->m_buffer = const_cast<T*>(data.data());
+ }
+
/**
* @brief Construct a new Array2D object from a Memory object
*
@@ -413,6 +443,21 @@ class Array3D : public Array1D<T> {
m_y_dim(y_dim),
m_x_dim(x_dim) {}
+ /**
+ * @brief Construct a new Array3D object from a Span
+ *
+ * @param data The Span that the Array will use.
+ * Must be avaible during the life time of the Array
+ * The user is responsible for allocation and deallocation.
+ */
+ Array3D(const aocommon::xt::Span<T, 3>& data)
+ : m_z_dim(data.shape(0)), m_y_dim(data.shape(1)), m_x_dim(data.shape(2)) {
+ this->m_memory = nullptr;
+ // This const_cast simplifies the transition to XTensor. It will
+ // be removed at the end of the transition, together with Array3D.
+ this->m_buffer = const_cast<T*>(data.data());
+ }
+
/**
* @brief Construct a new Array3D object from a Memory object
*
@@ -590,6 +635,24 @@ class Array4D : public Array1D<T> {
m_y_dim(y_dim),
m_x_dim(x_dim) {}
+ /**
+ * @brief Construct a new Array4D object from a Span
+ *
+ * @param data The Span that the Array will use.
+ * Must be avaible during the life time of the Array
+ * The user is responsible for allocation and deallocation.
+ */
+ Array4D(const aocommon::xt::Span<T, 4>& data)
+ : m_w_dim(data.shape(0)),
+ m_z_dim(data.shape(1)),
+ m_y_dim(data.shape(2)),
+ m_x_dim(data.shape(3)) {
+ this->m_memory = nullptr;
+ // This const_cast simplifies the transition to XTensor. It will
+ // be removed at the end of the transition, together with Array4D.
+ this->m_buffer = const_cast<T*>(data.data());
+ }
+
/**
* @brief Construct a new Array4D object from a Memory object
*
diff --git a/idg-lib/src/common/Plan.cpp b/idg-lib/src/common/Plan.cpp
index c7444876d2f41550cba53151b966b4d98c1aa6f1..171e9abdf6f4100354f6a0fcff24ed2bbef268a9 100644
--- a/idg-lib/src/common/Plan.cpp
+++ b/idg-lib/src/common/Plan.cpp
@@ -6,18 +6,26 @@
#include <algorithm> // max_element
#include <memory.h> // memcpy
+#include <xtensor/xview.hpp>
+#include <xtensor/xsort.hpp>
+
#include "Plan.h"
+#include "auxiliary.h"
using namespace std;
namespace idg {
Plan::Plan(const int kernel_size, const int subgrid_size, const int grid_size,
- const float cell_size, const Array1D<float>& shift,
- const Array1D<float>& frequencies, const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets, Options options)
- : m_subgrid_size(subgrid_size),
+ const float cell_size, const std::array<float, 2>& shift,
+ const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
+ Options options)
+ : m_shift(shift),
+ m_subgrid_size(subgrid_size),
m_cell_size(cell_size),
use_wtiles(false),
m_options(options) {
@@ -26,20 +34,21 @@ Plan::Plan(const int kernel_size, const int subgrid_size, const int grid_size,
#endif
WTiles dummy_wtiles;
- m_shift(0) = shift(0);
- m_shift(1) = shift(1);
initialize(kernel_size, subgrid_size, grid_size, cell_size, frequencies, uvw,
baselines, aterm_offsets, dummy_wtiles, options);
}
Plan::Plan(const int kernel_size, const int subgrid_size, const int grid_size,
- const float cell_size, const Array1D<float>& shift,
- const Array1D<float>& frequencies, const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets, WTiles& wtiles,
- Options options)
- : m_subgrid_size(subgrid_size),
+ const float cell_size, const std::array<float, 2>& shift,
+ const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
+ WTiles& wtiles, Options options)
+ : m_shift(shift),
+ m_subgrid_size(subgrid_size),
m_cell_size(cell_size),
use_wtiles(true),
m_options(options) {
@@ -47,8 +56,6 @@ Plan::Plan(const int kernel_size, const int subgrid_size, const int grid_size,
cout << "Plan::" << __func__ << " (with WTiles)" << endl;
#endif
- m_shift(0) = shift(0);
- m_shift(1) = shift(1);
initialize(kernel_size, subgrid_size, grid_size, cell_size, frequencies, uvw,
baselines, aterm_offsets, wtiles, options);
}
@@ -206,10 +213,11 @@ inline float meters_to_lambda(float meters, float frequency) {
std::vector<std::pair<int, int>> make_channel_groups(
float baseline_length, float uv_span_frequency, float image_size,
- const Array1D<float>& frequencies, unsigned int max_nr_channels = 0) {
+ const aocommon::xt::Span<float, 1>& frequencies,
+ unsigned int max_nr_channels = 0) {
std::vector<std::pair<int, int>> result;
- unsigned int nr_channels = frequencies.get_x_dim();
+ const size_t nr_channels = frequencies.size();
// There will be at most as many channel_groups as channels
result.reserve(nr_channels);
@@ -238,10 +246,11 @@ std::vector<std::pair<int, int>> make_channel_groups(
void Plan::initialize(
const int kernel_size, const int subgrid_size, const int grid_size,
- const float cell_size, const Array1D<float>& frequencies,
- const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets, WTiles& wtiles,
+ const float cell_size, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets, WTiles& wtiles,
const Options& options) {
#if defined(DEBUG)
cout << "Plan::" << __func__ << endl;
@@ -250,41 +259,33 @@ void Plan::initialize(
std::clog << "grid_size : " << grid_size << std::endl;
#endif
- // Check arguments
- assert(baselines.get_x_dim() == uvw.get_y_dim());
-
// Initialize arguments
- auto nr_baselines = uvw.get_y_dim();
- auto nr_timesteps = uvw.get_x_dim();
- auto nr_timeslots = aterm_offsets.get_x_dim() - 1;
- auto nr_channels = frequencies.get_x_dim();
- auto image_size = cell_size * grid_size; // TODO: remove
- auto wtile_size = wtiles.get_wtile_size();
+ const size_t nr_baselines = uvw.shape(0);
+ assert(baselines.size() == nr_baselines);
+ const size_t nr_timesteps = uvw.shape(1);
+ const size_t nr_timeslots = aterm_offsets.size() - 1;
+ const size_t nr_channels = frequencies.size();
+ const float image_size = cell_size * grid_size; // TODO: remove
+ const size_t wtile_size = wtiles.get_wtile_size();
// Get options
m_w_step = options.w_step;
- int nr_w_layers = options.nr_w_layers;
- int max_nr_timesteps_per_subgrid =
- min(options.max_nr_timesteps_per_subgrid, nr_timesteps);
- int max_nr_channels_per_subgrid = options.max_nr_channels_per_subgrid;
- bool plan_strict = options.plan_strict;
-
- // Spectral-line imaging
- bool simulate_spectral_line = options.simulate_spectral_line;
- auto nr_channels_ = nr_channels;
- if (simulate_spectral_line) {
- nr_channels = 1;
- }
+ const size_t nr_w_layers = options.nr_w_layers;
+ const size_t max_nr_timesteps_per_subgrid = min(
+ static_cast<size_t>(options.max_nr_timesteps_per_subgrid), nr_timesteps);
+ const size_t max_nr_channels_per_subgrid =
+ options.max_nr_channels_per_subgrid;
+ const bool plan_strict = options.plan_strict;
// Temporary metadata vector for individual baselines
- int max_nr_subgrids_per_baseline =
+ const size_t max_nr_subgrids_per_baseline =
max_nr_timesteps_per_subgrid > 0
? ((nr_timesteps + max_nr_timesteps_per_subgrid) /
max_nr_timesteps_per_subgrid) +
1
: nr_timesteps;
- idg::Array2D<Metadata> metadata_(nr_baselines,
- nr_channels * max_nr_subgrids_per_baseline);
+ xt::xtensor<Metadata, 2> metadata_(
+ {nr_baselines, nr_channels * max_nr_subgrids_per_baseline});
// Count the actual number of subgrids per baseline
std::vector<unsigned int> nr_subgrids_per_baseline(nr_baselines);
@@ -327,14 +328,14 @@ void Plan::initialize(
// Compute uv coordinates in pixels
struct DataPoint {
- unsigned timestep;
+ unsigned int timestep;
float u_pixels;
float v_pixels;
float w_lambda;
};
// Allocate datapoints for first and last channel in a group
- idg::Array2D<DataPoint> datapoints(nr_timesteps, 2);
+ xt::xtensor<DataPoint, 2> datapoints({nr_timesteps, 2});
for (auto channel_group : channel_groups) {
auto channel_begin = channel_group.first;
@@ -376,22 +377,22 @@ void Plan::initialize(
while (time_offset < nr_timesteps) {
// Create subgrid
subgrid.reset();
- int nr_timesteps_subgrid = 0;
+ size_t nr_timesteps_subgrid = 0;
// Load first visibility
- DataPoint first_datapoint = datapoints(time_offset, 0);
- const int first_timestep = first_datapoint.timestep;
+ const DataPoint first_datapoint = datapoints(time_offset, 0);
+ const size_t first_timestep = first_datapoint.timestep;
// Iterate all datapoints
for (; time_offset < nr_timesteps; time_offset++) {
// Visibility for first channel
- DataPoint visibility0 = datapoints(time_offset, 0);
+ const DataPoint visibility0 = datapoints(time_offset, 0);
const float u_pixels0 = visibility0.u_pixels;
const float v_pixels0 = visibility0.v_pixels;
const float w_lambda0 = visibility0.w_lambda;
// Visibility for last channel
- DataPoint visibility1 = datapoints(time_offset, 1);
+ const DataPoint visibility1 = datapoints(time_offset, 1);
const float u_pixels1 = visibility1.u_pixels;
const float v_pixels1 = visibility1.v_pixels;
@@ -431,7 +432,7 @@ void Plan::initialize(
}
// Compute time index for first visibility on subgrid
- auto time_index = bl * nr_timesteps + first_timestep;
+ const size_t time_index = bl * nr_timesteps + first_timestep;
// Finish subgrid
subgrid.finish();
@@ -439,8 +440,8 @@ void Plan::initialize(
// Add subgrid to metadata
if (subgrid.in_range()) {
Metadata m = {
- .time_index = (int)time_index, // time index
- .nr_timesteps = nr_timesteps_subgrid, // nr of timesteps
+ .time_index = static_cast<int>(time_index),
+ .nr_timesteps = static_cast<int>(nr_timesteps_subgrid),
.channel_begin = channel_begin,
.channel_end = channel_end,
.baseline = baseline, // baselines
@@ -454,18 +455,6 @@ void Plan::initialize(
unsigned subgrid_idx = nr_subgrids_per_baseline[bl];
metadata_(bl, subgrid_idx++) = m;
- // Add additional subgrids for subsequent frequencies
- if (simulate_spectral_line) {
- for (unsigned c = 1; c < nr_channels_; c++) {
- // Compute shifted subgrid for current frequency
- float shift = frequencies(c) / frequencies(0);
- Metadata m = metadata_(bl, subgrid_idx);
- m.coordinate.x *= shift;
- m.coordinate.y *= shift;
- metadata_(bl, subgrid_idx++) = m;
- }
- }
-
// Update number of subgrids
nr_subgrids_per_baseline[bl] = subgrid_idx;
} else if (plan_strict) {
@@ -530,7 +519,7 @@ void Plan::initialize(
} // end for bl
// Reserve aterm indices
- aterm_indices.resize(nr_baselines * nr_timesteps);
+ aterm_indices.resize({nr_baselines, nr_timesteps});
// Set aterm index for every timestep
#pragma omp parallel for
@@ -545,31 +534,31 @@ void Plan::initialize(
// The aterm index is equal to the timeslot
const unsigned int aterm_index = timeslot;
- // Determine number of timesteps in current aterm
- const unsigned nr_timesteps_per_aterm =
+ // Determine timesteps in current aterm
+ const unsigned int nr_timesteps_per_aterm =
next_aterm_offset - current_aterm_offset;
+ const unsigned int first_aterm = time_idx;
+ const unsigned int last_aterm = time_idx + nr_timesteps_per_aterm;
- for (unsigned timestep = 0; timestep < nr_timesteps_per_aterm;
- timestep++) {
- aterm_indices[bl * nr_timesteps + time_idx++] = aterm_index;
- }
+ xt::view(aterm_indices, bl, xt::range(first_aterm, last_aterm)) =
+ aterm_index;
+
+ time_idx += nr_timesteps_per_aterm;
}
}
- // Set nr_aterms
#pragma omp parallel for
- for (unsigned i = 0; i < metadata.size(); i++) {
- auto& m = metadata[i];
- auto aterm_index = aterm_indices[m.time_index];
- auto nr_aterms = 1;
- for (auto time = 0; time < m.nr_timesteps; time++) {
- auto aterm_index_current = aterm_indices[m.time_index + time];
+ for (size_t i = 0; i < metadata.size(); ++i) {
+ Metadata& m = metadata[i];
+ size_t aterm_index = aterm_indices[m.time_index];
+ size_t nr_aterms = 1;
+ for (size_t time = 0; time < static_cast<size_t>(m.nr_timesteps); ++time) {
+ const size_t aterm_index_current = aterm_indices[m.time_index + time];
if (aterm_index != aterm_index_current) {
- nr_aterms++;
+ ++nr_aterms;
aterm_index = aterm_index_current;
}
}
-
m.nr_aterms = nr_aterms;
}
@@ -678,8 +667,11 @@ void Plan::copy_metadata(void* ptr) const {
}
void Plan::copy_aterm_indices(void* ptr) const {
- memcpy(ptr, get_aterm_indices_ptr(),
- aterm_indices.size() * sizeof(unsigned int));
+ const size_t nr_baselines = aterm_indices.shape(0);
+ const size_t nr_timesteps = aterm_indices.shape(1);
+ const size_t sizeof_aterm_indices =
+ auxiliary::sizeof_aterm_indices(nr_baselines, nr_timesteps);
+ memcpy(ptr, aterm_indices.data(), sizeof_aterm_indices);
}
const unsigned int* Plan::get_aterm_indices_ptr(int bl) const {
diff --git a/idg-lib/src/common/Plan.h b/idg-lib/src/common/Plan.h
index 4846b5b97af21f128fc8f17b116a6f2a24362b93..3875b3950b54d4e8e3de3f3c8865f253e75f5942 100644
--- a/idg-lib/src/common/Plan.h
+++ b/idg-lib/src/common/Plan.h
@@ -10,10 +10,11 @@
#include <cmath>
#include <numeric>
#include <iterator>
+
+#include <aocommon/xt/span.h>
#include <omp.h>
#include "Types.h"
-#include "ArrayTypes.h"
#include "WTiles.h"
namespace idg {
@@ -53,10 +54,6 @@ class Plan {
// zero means no limit
unsigned max_nr_channels_per_subgrid = 0;
- // consider only first channel when creating subgrids,
- // add additional subgrids for every subsequent frequencies
- bool simulate_spectral_line = false;
-
// Imaging mode
Mode mode = Mode::FULL_POLARIZATION;
};
@@ -67,16 +64,21 @@ class Plan {
Plan(Plan&&) = default;
Plan(const int kernel_size, const int subgrid_size, const int grid_size,
- const float cell_size, const Array1D<float>& shift,
- const Array1D<float>& frequencies, const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets, Options options = Options());
+ const float cell_size, const std::array<float, 2>& shift,
+ const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
+ Options options = Options());
Plan(const int kernel_size, const int subgrid_size, const int grid_size,
- const float cell_size, const Array1D<float>& shift,
- const Array1D<float>& frequencies, const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets, WTiles& wtiles,
+ const float cell_size, const std::array<float, 2>& shift,
+ const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets, WTiles& wtiles,
Options options = Options());
// Destructor
@@ -84,10 +86,11 @@ class Plan {
void initialize(
const int kernel_size, const int subgrid_size, const int grid_size,
- const float cell_size, const Array1D<float>& frequencies,
- const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets, WTiles& wtiles,
+ const float cell_size, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets, WTiles& wtiles,
const Options& options);
// options
@@ -171,11 +174,11 @@ class Plan {
int get_subgrid_size() const { return m_subgrid_size; }
float get_w_step() const { return m_w_step; }
- const Array1D<float>& get_shift() const { return m_shift; }
+ const std::array<float, 2>& get_shift() const { return m_shift; }
float get_cell_size() const { return m_cell_size; }
private:
- Array1D<float> m_shift{2};
+ std::array<float, 2> m_shift;
int m_subgrid_size;
float m_w_step;
float m_cell_size;
@@ -185,7 +188,7 @@ class Plan {
std::vector<int> total_nr_visibilities_per_baseline;
WTileUpdateSet m_wtile_initialize_set;
WTileUpdateSet m_wtile_flush_set;
- std::vector<unsigned int> aterm_indices;
+ xt::xtensor<unsigned int, 2> aterm_indices;
bool use_wtiles;
Options m_options;
}; // class Plan
diff --git a/idg-lib/src/common/PlanC.cpp b/idg-lib/src/common/PlanC.cpp
index e56ad85fd3997ef98b4756df3540e23a51350aa2..e2b3a62e8f6c0deb7558a7128a60107bcd4be735 100644
--- a/idg-lib/src/common/PlanC.cpp
+++ b/idg-lib/src/common/PlanC.cpp
@@ -28,22 +28,28 @@ struct idg::Plan* Plan_init(int kernel_size, int subgrid_size, int grid_size,
int baselines_nr_baselines, int baselines_two,
unsigned int* aterm_offsets,
int aterm_offsets_nr_timeslots) {
- idg::Array1D<float> shift_array(3);
- shift_array(0) = 0;
- shift_array(1) = 0;
- shift_array(2) = 0;
- idg::Array1D<float> frequencies_array(frequencies, frequencies_nr_channels);
+ std::array<float, 2> shift{0.0f, 0.0f};
+ const std::array<size_t, 1> frequencies_shape{
+ static_cast<size_t>(frequencies_nr_channels)};
+ aocommon::xt::Span<float, 1> frequencies_span =
+ aocommon::xt::CreateSpan(frequencies, frequencies_shape);
assert(uvw_nr_coordinates == 3);
- idg::Array2D<idg::UVW<float>> uvw_array(uvw, uvw_nr_baselines,
- uvw_nr_timesteps);
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines_array(
- baselines, baselines_nr_baselines);
- idg::Array1D<unsigned int> aterm_offsets_array(aterm_offsets,
- aterm_offsets_nr_timeslots);
-
- return new idg::Plan(kernel_size, subgrid_size, grid_size, cell_size,
- shift_array, frequencies_array, uvw_array,
- baselines_array, aterm_offsets_array);
+ const std::array<size_t, 2> uvw_shape{static_cast<size_t>(uvw_nr_baselines),
+ static_cast<size_t>(uvw_nr_timesteps)};
+ aocommon::xt::Span<idg::UVW<float>, 2> uvw_span =
+ aocommon::xt::CreateSpan(uvw, uvw_shape);
+ const std::array<size_t, 1> baselines_shape{
+ static_cast<size_t>(baselines_nr_baselines)};
+ aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1> baselines_span =
+ aocommon::xt::CreateSpan(baselines, baselines_shape);
+ const std::array<size_t, 1> aterm_offsets_shape{
+ static_cast<size_t>(aterm_offsets_nr_timeslots)};
+ aocommon::xt::Span<unsigned int, 1> aterm_offsets_span =
+ aocommon::xt::CreateSpan(aterm_offsets, aterm_offsets_shape);
+
+ return new idg::Plan(kernel_size, subgrid_size, grid_size, cell_size, shift,
+ frequencies_span, uvw_span, baselines_span,
+ aterm_offsets_span);
}
} // extern "C"
diff --git a/idg-lib/src/common/Proxy.cpp b/idg-lib/src/common/Proxy.cpp
index 2be127b1f4ce2a573e88685aa11512be2fee68b8..8b0306f40eb3abc6018992747bbbb9fe66e36b4a 100644
--- a/idg-lib/src/common/Proxy.cpp
+++ b/idg-lib/src/common/Proxy.cpp
@@ -38,6 +38,27 @@ void Proxy::gridding(
aterm_offsets, spheroidal);
}
+void Proxy::gridding(
+ const Plan& plan, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<std::complex<float>, 4>& visibilities,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4>& aterms,
+ const aocommon::xt::Span<unsigned int, 1>& aterms_offsets,
+ const aocommon::xt::Span<float, 2>& taper) {
+ const Array1D<float> frequencies_array(frequencies);
+ const Array4D<std::complex<float>> visibilities_array(visibilities);
+ const Array2D<UVW<float>> uvw_array(uvw);
+ const Array1D<std::pair<unsigned int, unsigned int>> baselines_array(
+ baselines);
+ const Array4D<Matrix2x2<std::complex<float>>> aterms_array(aterms);
+ const Array1D<unsigned int> aterms_offsets_array(aterms_offsets);
+ const Array2D<float> taper_array(taper);
+ gridding(plan, frequencies_array, visibilities_array, uvw_array,
+ baselines_array, aterms_array, aterms_offsets_array, taper_array);
+}
+
void Proxy::degridding(
const Plan& plan, const Array1D<float>& frequencies,
Array4D<std::complex<float>>& visibilities, const Array2D<UVW<float>>& uvw,
@@ -60,6 +81,27 @@ void Proxy::degridding(
aterm_offsets, spheroidal);
}
+void Proxy::degridding(
+ const Plan& plan, const aocommon::xt::Span<float, 1>& frequencies,
+ aocommon::xt::Span<std::complex<float>, 4>& visibilities,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4>& aterms,
+ const aocommon::xt::Span<unsigned int, 1>& aterms_offsets,
+ const aocommon::xt::Span<float, 2>& taper) {
+ const Array1D<float> frequencies_array(frequencies);
+ Array4D<std::complex<float>> visibilities_array(visibilities);
+ const Array2D<UVW<float>> uvw_array(uvw);
+ const Array1D<std::pair<unsigned int, unsigned int>> baselines_array(
+ baselines);
+ const Array4D<Matrix2x2<std::complex<float>>> aterms_array(aterms);
+ const Array1D<unsigned int> aterms_offsets_array(aterms_offsets);
+ const Array2D<float> taper_array(taper);
+ degridding(plan, frequencies_array, visibilities_array, uvw_array,
+ baselines_array, aterms_array, aterms_offsets_array, taper_array);
+}
+
void Proxy::calibrate_init(
const unsigned int kernel_size, const Array2D<float>& frequencies,
Array4D<std::complex<float>>& visibilities, Array4D<float>& weights,
@@ -96,7 +138,7 @@ void Proxy::calibrate_init(
auto nr_timesteps = visibilities.get_z_dim();
auto nr_correlations = visibilities.get_x_dim();
assert(nr_correlations == 4);
- auto nr_baselines = baselines.get_x_dim();
+ const size_t nr_baselines = baselines.size();
auto nr_channel_blocks = frequencies.get_y_dim();
auto nr_channels_per_block = frequencies.get_x_dim();
@@ -188,14 +230,20 @@ void Proxy::calibrate_init(
plans[i].reserve(nr_channel_blocks);
for (unsigned int channel_block = 0; channel_block < nr_channel_blocks;
channel_block++) {
- Array1D<float> frequencies_channel_block(frequencies.data(channel_block),
- nr_channels_per_block);
+ const std::array<size_t, 1> frequencies_channel_block_shape{
+ nr_channels_per_block};
+ auto frequencies_channel_block = aocommon::xt::CreateSpan(
+ frequencies.data(channel_block), frequencies_channel_block_shape);
+ const std::array<size_t, 1> aterm_offsets_shape{aterm_offsets.size()};
+ auto aterm_offsets_span =
+ aocommon::xt::CreateSpan(aterm_offsets.data(), aterm_offsets_shape);
+ const std::array<size_t, 2> uvw_shape{nr_antennas - 1, nr_timesteps};
+ const std::array<size_t, 1> baselines_shape{nr_antennas - 1};
plans[i].push_back(make_plan(
kernel_size, frequencies_channel_block,
- Array2D<UVW<float>>(uvw1.data(i), nr_antennas - 1, nr_timesteps),
- Array1D<std::pair<unsigned int, unsigned int>>(baselines1.data(i),
- nr_antennas - 1),
- aterm_offsets, options));
+ aocommon::xt::CreateSpan(uvw1.data(i), uvw_shape),
+ aocommon::xt::CreateSpan(baselines1.data(i), baselines_shape),
+ aterm_offsets_span, options));
}
}
diff --git a/idg-lib/src/common/Proxy.h b/idg-lib/src/common/Proxy.h
index d07c088932bcbe5975b778380c616ac52352543c..48ae4779d522fb08a3c5e1ce45fe76872e83bb08 100644
--- a/idg-lib/src/common/Proxy.h
+++ b/idg-lib/src/common/Proxy.h
@@ -10,12 +10,16 @@
#include <cstring>
#include <utility> // pair
+#include <aocommon/xt/span.h>
+
+#include "ArrayTypes.h"
#include "RuntimeWrapper.h"
#include "ProxyInfo.h"
#include "Types.h"
#include "Plan.h"
#include "Report.h"
#include "Exception.h"
+#include "Tensor.h"
namespace idg {
enum DomainAtoDomainB {
@@ -81,6 +85,16 @@ class Proxy {
const Array1D<unsigned int>& aterm_offsets,
const Array2D<float>& taper);
+ void gridding(
+ const Plan& plan, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<std::complex<float>, 4>& visibilities,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4>& aterms,
+ const aocommon::xt::Span<unsigned int, 1>& aterms_offsets,
+ const aocommon::xt::Span<float, 2>& taper);
+
/**
* @brief Degrid (predict) visibilities, applying A-terms.
*
@@ -120,6 +134,16 @@ class Proxy {
const Array4D<Matrix2x2<std::complex<float>>>& aterms,
const Array1D<unsigned int>& aterm_offsets, const Array2D<float>& taper);
+ void degridding(
+ const Plan& plan, const aocommon::xt::Span<float, 1>& frequencies,
+ aocommon::xt::Span<std::complex<float>, 4>& visibilities,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4>& aterms,
+ const aocommon::xt::Span<unsigned int, 1>& aterms_offsets,
+ const aocommon::xt::Span<float, 2>& taper);
+
/**
* @brief Prepare a calibration cycle
*
@@ -276,6 +300,37 @@ class Proxy {
return Array4D<T>(allocate_memory(bytes), d_dim, c_dim, b_dim, a_dim);
};
+ template <typename T, size_t Dimensions>
+ Tensor<T, Dimensions> allocate_tensor(
+ const std::initializer_list<size_t> shape) {
+ assert(shape.size() == Dimensions);
+ std::array<size_t, Dimensions> shape_array;
+ size_t bytes = sizeof(T);
+ for (size_t i = 0; i < shape.size(); i++) {
+ const size_t dimension = *(shape.begin() + i);
+ shape_array[i] = dimension;
+ bytes *= dimension;
+ }
+ return Tensor<T, Dimensions>(allocate_memory(bytes), shape_array);
+ }
+
+ template <typename T, size_t Dimensions>
+ aocommon::xt::Span<T, Dimensions> allocate_span(
+ const std::initializer_list<size_t> shape) {
+ assert(shape.size() == Dimensions);
+ std::array<size_t, Dimensions> shape_array;
+ size_t bytes = sizeof(T);
+ for (size_t i = 0; i < shape.size(); i++) {
+ const size_t dimension = *(shape.begin() + i);
+ shape_array[i] = dimension;
+ bytes *= dimension;
+ }
+ std::unique_ptr<auxiliary::Memory> memory = allocate_memory(bytes);
+ T* ptr = reinterpret_cast<T*>(memory->data());
+ memory_.push_back(std::move(memory));
+ return aocommon::xt::CreateSpan(ptr, shape_array);
+ }
+
//! Methods for grid management
virtual std::shared_ptr<Grid> allocate_grid(size_t nr_w_layers,
size_t nr_correlations,
@@ -320,12 +375,11 @@ class Proxy {
* @param shift
*/
virtual void init_cache(int subgrid_size, float cell_size, float w_step,
- const Array1D<float>& shift) {
+ const std::array<float, 2>& shift) {
m_cache_state.subgrid_size = subgrid_size;
m_cache_state.cell_size = cell_size;
m_cache_state.w_step = w_step;
- m_cache_state.shift(0) = shift(0);
- m_cache_state.shift(1) = shift(1);
+ m_cache_state.shift = shift;
};
// The cache needs to have been initialized by call to init_cache first
@@ -351,10 +405,11 @@ class Proxy {
* @return std::unique_ptr<Plan>
*/
virtual std::unique_ptr<Plan> make_plan(
- const int kernel_size, const Array1D<float>& frequencies,
- const Array2D<UVW<float>>& uvw,
- const Array1D<std::pair<unsigned int, unsigned int>>& baselines,
- const Array1D<unsigned int>& aterm_offsets,
+ const int kernel_size, const aocommon::xt::Span<float, 1>& frequencies,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>&
+ baselines,
+ const aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
Plan::Options options = Plan::Options()) {
options.w_step = m_cache_state.w_step;
return std::unique_ptr<Plan>(
@@ -467,9 +522,14 @@ class Proxy {
int subgrid_size;
float cell_size;
float w_step;
- Array1D<float> shift{2};
+ std::array<float, 2> shift;
} m_cache_state;
+ void free_memory() { memory_.clear(); };
+
+ private:
+ std::vector<std::unique_ptr<auxiliary::Memory>> memory_;
+
}; // end class Proxy
} // namespace proxy
diff --git a/idg-lib/src/common/ProxyC.cpp b/idg-lib/src/common/ProxyC.cpp
index 7d632db7330b711643060b3fcfa64484b9066a55..fb6010dd7e876e78f71208d85a29a2a834d6e63b 100644
--- a/idg-lib/src/common/ProxyC.cpp
+++ b/idg-lib/src/common/ProxyC.cpp
@@ -44,17 +44,35 @@ void Proxy_gridding(struct Proxy* p, int kernel_size, int subgrid_size,
std::pair<unsigned int, unsigned int>* baselines,
std::complex<float>* aterms, unsigned int* aterm_offsets,
float* taper) {
- idg::Array1D<float> frequencies_(frequencies, nr_channels);
- idg::Array4D<std::complex<float>> visibilities_(
- visibilities, nr_baselines, nr_timesteps, nr_channels, nr_correlations);
- idg::Array2D<idg::UVW<float>> uvw_(uvw, nr_baselines, nr_timesteps);
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines_(baselines,
- nr_baselines);
- idg::Array4D<idg::Matrix2x2<std::complex<float>>> aterms_(
- (idg::Matrix2x2<std::complex<float>>*)aterms, nr_timeslots, nr_stations,
- subgrid_size, subgrid_size);
- idg::Array1D<unsigned int> aterm_offsets_(aterm_offsets, nr_timeslots + 1);
- idg::Array2D<float> taper_(taper, subgrid_size, subgrid_size);
+ const std::array<size_t, 1> frequencies_shape{
+ static_cast<size_t>(nr_channels)};
+ const std::array<size_t, 2> uvw_shape{static_cast<size_t>(nr_baselines),
+ static_cast<size_t>(nr_timesteps)};
+ const std::array<size_t, 4> visibilities_shape{
+ static_cast<size_t>(nr_baselines), static_cast<size_t>(nr_timesteps),
+ static_cast<size_t>(nr_channels), static_cast<size_t>(nr_correlations)};
+ const std::array<size_t, 1> baselines_shape{
+ static_cast<size_t>(nr_baselines)};
+ const std::array<size_t, 4> aterms_shape{
+ static_cast<size_t>(nr_timeslots), static_cast<size_t>(nr_stations),
+ static_cast<size_t>(subgrid_size), static_cast<size_t>(subgrid_size)};
+ const std::array<size_t, 1> aterm_offsets_shape{
+ static_cast<size_t>(nr_timeslots + 1)};
+ const std::array<size_t, 2> taper_shape{static_cast<size_t>(subgrid_size),
+ static_cast<size_t>(subgrid_size)};
+
+ auto frequencies_span =
+ aocommon::xt::CreateSpan(frequencies, frequencies_shape);
+ auto uvw_span = aocommon::xt::CreateSpan(uvw, uvw_shape);
+ auto visibilities_span =
+ aocommon::xt::CreateSpan(visibilities, visibilities_shape);
+ auto baselines_span = aocommon::xt::CreateSpan(baselines, baselines_shape);
+ auto aterms_span = aocommon::xt::CreateSpan(
+ reinterpret_cast<idg::Matrix2x2<std::complex<float>>*>(aterms),
+ aterms_shape);
+ auto aterm_offsets_span =
+ aocommon::xt::CreateSpan(aterm_offsets, aterm_offsets_shape);
+ auto taper_span = aocommon::xt::CreateSpan(taper, taper_shape);
idg::Plan::Options options;
options.mode = nr_correlations == 4 ? idg::Plan::Mode::FULL_POLARIZATION
@@ -62,12 +80,12 @@ void Proxy_gridding(struct Proxy* p, int kernel_size, int subgrid_size,
std::unique_ptr<idg::Plan> plan = ExitOnException(
&idg::proxy::Proxy::make_plan, reinterpret_cast<idg::proxy::Proxy*>(p),
- kernel_size, frequencies_, uvw_, baselines_, aterm_offsets_, options);
+ kernel_size, frequencies_span, uvw_span, baselines_span,
+ aterm_offsets_span, options);
- ExitOnException(&idg::proxy::Proxy::gridding,
- reinterpret_cast<idg::proxy::Proxy*>(p), *plan, frequencies_,
- visibilities_, uvw_, baselines_, aterms_, aterm_offsets_,
- taper_);
+ reinterpret_cast<idg::proxy::Proxy*>(p)->gridding(
+ *plan, frequencies_span, visibilities_span, uvw_span, baselines_span,
+ aterms_span, aterm_offsets_span, taper_span);
}
void Proxy_degridding(struct Proxy* p, int kernel_size, int subgrid_size,
@@ -78,17 +96,35 @@ void Proxy_degridding(struct Proxy* p, int kernel_size, int subgrid_size,
std::pair<unsigned int, unsigned int>* baselines,
std::complex<float>* aterms, unsigned int* aterm_offsets,
float* taper) {
- idg::Array1D<float> frequencies_(frequencies, nr_channels);
- idg::Array4D<std::complex<float>> visibilities_(
- visibilities, nr_baselines, nr_timesteps, nr_channels, nr_correlations);
- idg::Array2D<idg::UVW<float>> uvw_(uvw, nr_baselines, nr_timesteps);
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines_(baselines,
- nr_baselines);
- idg::Array4D<idg::Matrix2x2<std::complex<float>>> aterms_(
- (idg::Matrix2x2<std::complex<float>>*)aterms, nr_timeslots, nr_stations,
- subgrid_size, subgrid_size);
- idg::Array1D<unsigned int> aterm_offsets_(aterm_offsets, nr_timeslots + 1);
- idg::Array2D<float> taper_(taper, subgrid_size, subgrid_size);
+ const std::array<size_t, 1> frequencies_shape{
+ static_cast<size_t>(nr_channels)};
+ const std::array<size_t, 2> uvw_shape{static_cast<size_t>(nr_baselines),
+ static_cast<size_t>(nr_timesteps)};
+ const std::array<size_t, 4> visibilities_shape{
+ static_cast<size_t>(nr_baselines), static_cast<size_t>(nr_timesteps),
+ static_cast<size_t>(nr_channels), static_cast<size_t>(nr_correlations)};
+ const std::array<size_t, 1> baselines_shape{
+ static_cast<size_t>(nr_baselines)};
+ const std::array<size_t, 4> aterms_shape{
+ static_cast<size_t>(nr_timeslots), static_cast<size_t>(nr_stations),
+ static_cast<size_t>(subgrid_size), static_cast<size_t>(subgrid_size)};
+ const std::array<size_t, 1> aterm_offsets_shape{
+ static_cast<size_t>(nr_timeslots + 1)};
+ const std::array<size_t, 2> taper_shape{static_cast<size_t>(subgrid_size),
+ static_cast<size_t>(subgrid_size)};
+
+ auto frequencies_span =
+ aocommon::xt::CreateSpan(frequencies, frequencies_shape);
+ auto uvw_span = aocommon::xt::CreateSpan(uvw, uvw_shape);
+ auto visibilities_span =
+ aocommon::xt::CreateSpan(visibilities, visibilities_shape);
+ auto baselines_span = aocommon::xt::CreateSpan(baselines, baselines_shape);
+ auto aterms_span = aocommon::xt::CreateSpan(
+ reinterpret_cast<idg::Matrix2x2<std::complex<float>>*>(aterms),
+ aterms_shape);
+ auto aterm_offsets_span =
+ aocommon::xt::CreateSpan(aterm_offsets, aterm_offsets_shape);
+ auto taper_span = aocommon::xt::CreateSpan(taper, taper_shape);
idg::Plan::Options options;
options.mode = nr_correlations == 4 ? idg::Plan::Mode::FULL_POLARIZATION
@@ -96,20 +132,20 @@ void Proxy_degridding(struct Proxy* p, int kernel_size, int subgrid_size,
std::unique_ptr<idg::Plan> plan = ExitOnException(
&idg::proxy::Proxy::make_plan, reinterpret_cast<idg::proxy::Proxy*>(p),
- kernel_size, frequencies_, uvw_, baselines_, aterm_offsets_, options);
+ kernel_size, frequencies_span, uvw_span, baselines_span,
+ aterm_offsets_span, options);
- ExitOnException(&idg::proxy::Proxy::degridding,
- reinterpret_cast<idg::proxy::Proxy*>(p), *plan, frequencies_,
- visibilities_, uvw_, baselines_, aterms_, aterm_offsets_,
- taper_);
+ reinterpret_cast<idg::proxy::Proxy*>(p)->degridding(
+ *plan, frequencies_span, visibilities_span, uvw_span, baselines_span,
+ aterms_span, aterm_offsets_span, taper_span);
}
void Proxy_init_cache(struct Proxy* p, unsigned int subgrid_size,
const float cell_size, float w_step, float* shift) {
- idg::Array1D<float> shift_(shift, 2);
+ const std::array<float, 2> shift_array{shift[0], shift[1]};
ExitOnException(&idg::proxy::Proxy::init_cache,
reinterpret_cast<idg::proxy::Proxy*>(p), subgrid_size,
- cell_size, w_step, shift_);
+ cell_size, w_step, shift_array);
}
void Proxy_calibrate_init(struct Proxy* p, unsigned int kernel_size,
@@ -131,8 +167,11 @@ void Proxy_calibrate_init(struct Proxy* p, unsigned int kernel_size,
nr_timesteps, nr_channels, nr_correlations);
idg::Array2D<idg::UVW<float>> uvw_((idg::UVW<float>*)uvw, nr_baselines,
nr_timesteps);
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines_(
- (std::pair<unsigned int, unsigned int>*)baselines, nr_baselines);
+ const std::array<size_t, 1> baselines_shape{
+ static_cast<size_t>(nr_baselines)};
+ auto baselines_ = aocommon::xt::CreateSpan(
+ reinterpret_cast<std::pair<unsigned int, unsigned int>*>(baselines),
+ baselines_shape);
idg::Array1D<unsigned int> aterm_offsets_(aterm_offsets, nr_timeslots + 1);
idg::Array2D<float> spheroidal_(spheroidal, subgrid_size, subgrid_size);
diff --git a/idg-lib/src/common/Tensor.h b/idg-lib/src/common/Tensor.h
new file mode 100644
index 0000000000000000000000000000000000000000..bcf6c5d8f0068b6388c178b8f7cbb477a7d79be3
--- /dev/null
+++ b/idg-lib/src/common/Tensor.h
@@ -0,0 +1,27 @@
+// Copyright (C) 2023 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include <aocommon/xt/span.h>
+
+#include "auxiliary.h"
+
+namespace idg {
+
+/// Structure that holds a span and the memory allocated for that span.
+template <typename T, size_t Dimensions>
+class Tensor {
+ public:
+ Tensor(std::unique_ptr<auxiliary::Memory> memory,
+ const std::array<size_t, Dimensions>& shape)
+ : memory_(std::move(memory)),
+ span_(aocommon::xt::CreateSpan<T, Dimensions>(
+ reinterpret_cast<T*>(memory_->data()), shape)){};
+
+ aocommon::xt::Span<T, Dimensions>& Span() { return span_; }
+
+ private:
+ std::unique_ptr<auxiliary::Memory> memory_;
+ aocommon::xt::Span<T, Dimensions> span_;
+};
+
+} // namespace idg
\ No newline at end of file
diff --git a/idg-lib/src/common/WTiles.cpp b/idg-lib/src/common/WTiles.cpp
index fd5a20b995b9ed74d9c1d3354218f91aa2b799ff..287c2422b0adaa7f4b9c0d91048a4f55eef24dc4 100644
--- a/idg-lib/src/common/WTiles.cpp
+++ b/idg-lib/src/common/WTiles.cpp
@@ -116,7 +116,7 @@ std::vector<int> compute_w_padded_tile_sizes(const idg::Coordinate* coordinates,
int compute_w_padded_tile_size_max(const WTileUpdateSet& wtile_set,
const int tile_size, const int subgrid_size,
const float image_size, const float w_step,
- const float shift_l, const float shift_m) {
+ const std::array<float, 2>& shift) {
int w_padded_tile_size_max = 0;
for (unsigned int i = 0; i < wtile_set.size(); i++) {
@@ -126,7 +126,7 @@ int compute_w_padded_tile_size_max(const WTileUpdateSet& wtile_set,
wtile_info.wtile_coordinates;
const int padded_tile_size = tile_size + subgrid_size;
const float image_size_shift =
- image_size + 2 * std::max(std::abs(shift_l), std::abs(shift_m));
+ image_size + 2 * std::max(std::abs(shift[0]), std::abs(shift[1]));
const std::vector<int> w_padded_tile_sizes = compute_w_padded_tile_sizes(
tile_coordinates.data(), nr_tiles, w_step, image_size, image_size_shift,
padded_tile_size);
diff --git a/idg-lib/src/common/WTiles.h b/idg-lib/src/common/WTiles.h
index c9df5c11e73073cc047907db551e8f84bcb369b1..ab921977cad002edcce8ede6e3a8eea495e834dd 100644
--- a/idg-lib/src/common/WTiles.h
+++ b/idg-lib/src/common/WTiles.h
@@ -191,7 +191,7 @@ std::vector<int> compute_w_padded_tile_sizes(const idg::Coordinate* coordinates,
int compute_w_padded_tile_size_max(const WTileUpdateSet& wtile_set,
const int tile_size, const int subgrid_size,
const float image_size, const float w_step,
- const float shift_l, const float shift_m);
+ const std::array<float, 2>& shift);
} // namespace idg
#endif
diff --git a/idg-util/src/initialize/Data.cpp b/idg-util/src/initialize/Data.cpp
index 0c96fadf37d20e880e0030b8b5a6d46f01373cd0..ca6688d4d0fed1ca62984418396604404d25efc8 100644
--- a/idg-util/src/initialize/Data.cpp
+++ b/idg-util/src/initialize/Data.cpp
@@ -215,40 +215,33 @@ void Data::limit_nr_baselines(unsigned int n) {
std::swap(m_baselines, baselines_selected);
}
-void Data::get_frequencies(Array1D<float>& frequencies, float image_size,
+void Data::get_frequencies(aocommon::xt::Span<float, 1>& frequencies,
+ float image_size,
unsigned int channel_offset) const {
- auto nr_channels = frequencies.get_x_dim();
- auto max_uv = get_max_uv();
- float frequency_increment = SPEED_OF_LIGHT / (max_uv * image_size);
- for (unsigned chan = 0; chan < nr_channels; chan++) {
+ const size_t nr_channels = frequencies.size();
+ const float max_uv = get_max_uv();
+ const float frequency_increment = SPEED_OF_LIGHT / (max_uv * image_size);
+ for (size_t chan = 0; chan < nr_channels; chan++) {
frequencies(chan) =
start_frequency + frequency_increment * (chan + channel_offset);
}
}
-Array2D<UVW<float>> Data::get_uvw(unsigned int nr_baselines,
- unsigned int nr_timesteps,
- unsigned int baseline_offset,
- unsigned int time_offset,
- float integration_time) const {
- Array2D<UVW<float>> uvw(nr_baselines, nr_timesteps);
- get_uvw(uvw, baseline_offset, time_offset, integration_time);
- return uvw;
-}
-
-void Data::get_uvw(Array2D<UVW<float>>& uvw, unsigned int baseline_offset,
- unsigned int time_offset, float integration_time) const {
- unsigned int nr_baselines_total = m_baselines.size();
- unsigned int nr_baselines = uvw.get_y_dim();
- unsigned int nr_timesteps = uvw.get_x_dim();
+void Data::get_uvw(aocommon::xt::Span<UVW<float>, 2>& uvw,
+ unsigned int baseline_offset, unsigned int time_offset,
+ float integration_time) const {
+ const size_t nr_baselines_total = m_baselines.size();
+ const size_t nr_baselines = uvw.shape(0);
+ const size_t nr_timesteps = uvw.shape(1);
if (baseline_offset + nr_baselines > nr_baselines_total) {
- std::cerr << "Out-of-range baselines selected: ";
+ std::stringstream message;
+ message << "Out-of-range baselines selected: ";
if (baseline_offset > 0) {
- std::cerr << baseline_offset << " + " << nr_baselines;
+ message << baseline_offset << " + " << nr_baselines;
}
- std::cerr << nr_baselines << " > " << nr_baselines_total << std::endl;
- nr_baselines = nr_baselines_total;
+ message << nr_baselines << " > " << nr_baselines_total << std::endl;
+ throw std::runtime_error(message.str());
}
// Evaluate uvw per baseline
diff --git a/idg-util/src/initialize/Data.h b/idg-util/src/initialize/Data.h
index 147ed61a444239df803465fb31ad7b7c6411d176..4183f4f2d6605297808478c05e2f9e4dac8ff5ca 100644
--- a/idg-util/src/initialize/Data.h
+++ b/idg-util/src/initialize/Data.h
@@ -47,17 +47,11 @@ class Data {
unsigned int get_nr_baselines() const { return m_baselines.size(); };
- void get_frequencies(Array1D<float>& frequencies, float image_size,
- unsigned int channel_offset = 0) const;
+ void get_frequencies(aocommon::xt::Span<float, 1>& frequencies,
+ float image_size, unsigned int channel_offset = 0) const;
- Array2D<UVW<float>> get_uvw(unsigned int nr_baselines,
- unsigned int nr_timesteps,
- unsigned int baseline_offset = 0,
- unsigned int time_offset = 0,
- float integration_time = 1.0f) const;
-
- void get_uvw(Array2D<UVW<float>>& uvw, unsigned int baseline_offset = 0,
- unsigned int time_offset = 0,
+ void get_uvw(aocommon::xt::Span<UVW<float>, 2>& uvw,
+ unsigned int baseline_offset = 0, unsigned int time_offset = 0,
float integration_time = 1.0f) const;
float get_max_uv() const;
diff --git a/idg-util/src/initialize/datac.h b/idg-util/src/initialize/datac.h
index c7bb3bef79483583910dc6507d86c4a4f979e9b3..6381c317318099d41e30ad91a267cb0d59b4bbe0 100644
--- a/idg-util/src/initialize/datac.h
+++ b/idg-util/src/initialize/datac.h
@@ -41,15 +41,18 @@ float DATA_get_nr_baselines(idg::Data* data) {
void DATA_get_frequencies(idg::Data* data, void* ptr, unsigned int nr_channels,
float image_size, unsigned int channel_offset) {
- idg::Array1D<float> frequencies((float*)ptr, nr_channels);
+ const std::array<size_t, 1> frequencies_shape{nr_channels};
+ auto frequencies = aocommon::xt::CreateSpan(reinterpret_cast<float*>(ptr),
+ frequencies_shape);
data->get_frequencies(frequencies, image_size, channel_offset);
}
void DATA_get_uvw(idg::Data* data, void* ptr, unsigned int nr_baselines,
unsigned int nr_timesteps, unsigned int baseline_offset,
unsigned int time_offset, float integration_time) {
- idg::Array2D<idg::UVW<float>> uvw((idg::UVW<float>*)ptr, nr_baselines,
- nr_timesteps);
+ const std::array<size_t, 2> uvw_shape{nr_baselines, nr_timesteps};
+ auto uvw = aocommon::xt::CreateSpan(reinterpret_cast<idg::UVW<float>*>(ptr),
+ uvw_shape);
data->get_uvw(uvw, baseline_offset, time_offset, integration_time);
}
diff --git a/idg-util/src/initialize/init.cpp b/idg-util/src/initialize/init.cpp
index 83a4e713c1f181e8bafd450d66187cecbb0d106a..ee169f3e4e1c7d186f7535bd4aee81173325b368 100644
--- a/idg-util/src/initialize/init.cpp
+++ b/idg-util/src/initialize/init.cpp
@@ -3,8 +3,8 @@
#include "init.h"
-namespace idg {
-
+#include <xtensor/xview.hpp>
+namespace {
// Function to compute spheroidal. Based on reference code by BvdT.
float evaluate_spheroidal(float nu) {
float P[2][5] = {
@@ -47,341 +47,176 @@ float evaluate_spheroidal(float nu) {
return (1.0 - nusq) * (top / bot);
}
}
-
-// TODO: make generic, not spheroidal specific
-// TODO: use real-to-complex and complex-to-real FFT
-void resize_spheroidal(float* __restrict__ spheroidal_in, unsigned int size_in,
- float* __restrict__ spheroidal_out,
- unsigned int size_out) {
- auto in_ft = new std::complex<float>[size_in * size_in];
- auto out_ft = new std::complex<float>[size_out * size_out];
-
- for (unsigned i = 0; i < size_in; i++) {
- for (unsigned j = 0; j < size_in; j++) {
- in_ft[i * size_in + j] = spheroidal_in[i * size_in + j];
- }
- }
- idg::fft2f(size_in, in_ft);
-
- int offset = int((size_out - size_in) / 2);
-
- for (unsigned i = 0; i < size_in; i++) {
- for (unsigned j = 0; j < size_in; j++) {
- out_ft[(i + offset) * size_out + (j + offset)] = in_ft[i * size_in + j];
- }
- }
- idg::ifft2f(size_out, out_ft);
-
- float s = 1.0f / (size_in * size_in);
- for (unsigned i = 0; i < size_out; i++) {
- for (unsigned j = 0; j < size_out; j++) {
- spheroidal_out[i * size_out + j] = out_ft[i * size_out + j].real() * s;
- }
- }
-
- delete[] in_ft;
- delete[] out_ft;
-}
-
-void init_example_aterms(Array4D<Matrix2x2<std::complex<float>>>& aterms) {
- unsigned int nr_timeslots = aterms.get_w_dim();
- unsigned int nr_stations = aterms.get_z_dim();
- unsigned int height = aterms.get_y_dim();
- unsigned int width = aterms.get_x_dim();
-
- for (unsigned t = 0; t < nr_timeslots; t++) {
- for (unsigned ant = 0; ant < nr_stations; ant++) {
- for (unsigned y = 0; y < height; y++) {
- for (unsigned x = 0; x < width; x++) {
- float scale = ((float)(t + 1) / nr_timeslots);
- std::complex<float> valueXX = std::complex<float>(scale * 1.0, 1.1);
- std::complex<float> valueXY = std::complex<float>(scale * 0.8, 0.9);
- std::complex<float> valueYX = std::complex<float>(scale * 0.6, 1.7);
- std::complex<float> valueYY = std::complex<float>(scale * 0.4, 0.5);
- const Matrix2x2<std::complex<float>> aterm = {valueXX, valueXY,
- valueYX, valueYY};
- aterms(t, ant, y, x) = aterm;
- }
- }
- }
- }
-}
+} // namespace
+namespace idg {
/*
* Memory-allocation is handled by Proxy
*/
-Array1D<float> get_example_frequencies(proxy::Proxy& proxy,
- unsigned int nr_channels,
- float start_frequency,
- float frequency_increment) {
- using T = float;
- Array1D<T> frequencies = proxy.allocate_array1d<T>(nr_channels);
-
- for (unsigned chan = 0; chan < nr_channels; chan++) {
- frequencies(chan) = start_frequency + frequency_increment * chan;
- }
-
+aocommon::xt::Span<float, 1> get_example_frequencies(
+ proxy::Proxy& proxy, unsigned int nr_channels, float start_frequency,
+ float frequency_increment) {
+ auto frequencies = proxy.allocate_span<float, 1>({nr_channels});
+ init_example_frequencies(frequencies, start_frequency, frequency_increment);
return frequencies;
}
-Array4D<std::complex<float>> get_dummy_visibilities(
+aocommon::xt::Span<std::complex<float>, 4> get_dummy_visibilities(
proxy::Proxy& proxy, unsigned int nr_baselines, unsigned int nr_timesteps,
unsigned int nr_channels, unsigned int nr_correlations) {
assert(nr_correlations == 2 || nr_correlations == 4);
-
- Array4D<std::complex<float>> visibilities =
- proxy.allocate_array4d<std::complex<float>>(nr_baselines, nr_timesteps,
- nr_channels, nr_correlations);
-
-// Set all visibilities
-#pragma omp parallel for
- for (unsigned bl = 0; bl < nr_baselines; bl++) {
- for (unsigned time = 0; time < nr_timesteps; time++) {
- for (unsigned chan = 0; chan < nr_channels; chan++) {
- if (nr_correlations == 2) {
- visibilities(bl, time, chan, 0) = 1.0f;
- visibilities(bl, time, chan, 1) = 1.0f;
- } else if (nr_correlations == 4) {
- visibilities(bl, time, chan, 0) = 1.0f;
- visibilities(bl, time, chan, 1) = 0.0f;
- visibilities(bl, time, chan, 2) = 0.0f;
- visibilities(bl, time, chan, 3) = 1.0f;
- }
- }
- }
- }
-
+ auto visibilities = proxy.allocate_span<std::complex<float>, 4>(
+ {nr_baselines, nr_timesteps, nr_channels, nr_correlations});
+ init_dummy_visibilities(visibilities);
return visibilities;
}
-Array4D<std::complex<float>> get_example_visibilities(
- proxy::Proxy& proxy, Array2D<UVW<float>>& uvw, Array1D<float>& frequencies,
- float image_size, unsigned int nr_correlations, unsigned int grid_size,
- unsigned int nr_point_sources, int max_pixel_offset,
+aocommon::xt::Span<std::complex<float>, 4> get_example_visibilities(
+ proxy::Proxy& proxy, aocommon::xt::Span<UVW<float>, 2>& uvw,
+ aocommon::xt::Span<float, 1>& frequencies, float image_size,
+ unsigned int nr_correlations, unsigned int grid_size,
+ unsigned int nr_point_sources, unsigned int max_pixel_offset,
unsigned int random_seed, float amplitude) {
- unsigned int nr_baselines = uvw.get_y_dim();
- unsigned int nr_timesteps = uvw.get_x_dim();
- unsigned int nr_channels = frequencies.get_x_dim();
-
- Array4D<std::complex<float>> visibilities =
- proxy.allocate_array4d<std::complex<float>>(nr_baselines, nr_timesteps,
- nr_channels, nr_correlations);
-
- if (max_pixel_offset == -1) {
- max_pixel_offset = grid_size / 2;
- }
-
- srand(random_seed);
-
- for (unsigned i = 0; i < nr_point_sources; i++) {
- float x_offset = static_cast<float>(random()) / RAND_MAX;
- float y_offset = static_cast<float>(random()) / RAND_MAX;
- x_offset = (x_offset * (max_pixel_offset)) - (max_pixel_offset / 2);
- y_offset = (y_offset * (max_pixel_offset)) - (max_pixel_offset / 2);
-
- add_pt_src(visibilities, uvw, frequencies, image_size, grid_size, x_offset,
- y_offset, amplitude);
- }
-
+ const size_t nr_baselines = uvw.shape(0);
+ const size_t nr_timesteps = uvw.shape(1);
+ const size_t nr_channels = frequencies.shape(0);
+ auto visibilities = proxy.allocate_span<std::complex<float>, 4>(
+ {nr_baselines, nr_timesteps, nr_channels, nr_correlations});
+ init_example_visibilities(visibilities, uvw, frequencies, image_size,
+ grid_size, nr_point_sources, max_pixel_offset,
+ random_seed, amplitude);
return visibilities;
}
-Array1D<std::pair<unsigned int, unsigned int>> get_example_baselines(
- proxy::Proxy& proxy, unsigned int nr_stations, unsigned int nr_baselines) {
+aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>
+get_example_baselines(proxy::Proxy& proxy, unsigned int nr_stations,
+ unsigned int nr_baselines) {
using T = std::pair<unsigned int, unsigned int>;
- Array1D<T> baselines = proxy.allocate_array1d<T>(nr_baselines);
-
- unsigned bl = 0;
-
- for (unsigned station1 = 0; station1 < nr_stations; station1++) {
- for (unsigned station2 = station1 + 1; station2 < nr_stations; station2++) {
- if (bl >= nr_baselines) {
- break;
- }
- baselines(bl) = std::pair<unsigned int, unsigned int>(station1, station2);
- bl++;
- }
- }
-
+ auto baselines = proxy.allocate_span<T, 1>({nr_baselines});
+ init_example_baselines(baselines, nr_stations);
return baselines;
}
-Array4D<Matrix2x2<std::complex<float>>> get_identity_aterms(
+aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4> get_identity_aterms(
proxy::Proxy& proxy, unsigned int nr_timeslots, unsigned int nr_stations,
unsigned int height, unsigned int width) {
using T = Matrix2x2<std::complex<float>>;
- Array4D<T> aterms =
- proxy.allocate_array4d<T>(nr_timeslots, nr_stations, height, width);
-
- for (unsigned t = 0; t < nr_timeslots; t++) {
- for (unsigned ant = 0; ant < nr_stations; ant++) {
- for (unsigned y = 0; y < height; y++) {
- for (unsigned x = 0; x < width; x++) {
- std::complex<float> valueXX = std::complex<float>(1.0, 0.0);
- std::complex<float> valueXY = std::complex<float>(0.0, 0.0);
- std::complex<float> valueYX = std::complex<float>(0.0, 0.0);
- std::complex<float> valueYY = std::complex<float>(1.0, 0.0);
- const Matrix2x2<std::complex<float>> aterm = {valueXX, valueXY,
- valueYX, valueYY};
- aterms(t, ant, y, x) = aterm;
- }
- }
- }
- }
-
+ auto aterms =
+ proxy.allocate_span<T, 4>({nr_timeslots, nr_stations, height, width});
+ init_identity_aterms(aterms);
return aterms;
}
-Array4D<Matrix2x2<std::complex<float>>> get_example_aterms(
+aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4> get_example_aterms(
proxy::Proxy& proxy, unsigned int nr_timeslots, unsigned int nr_stations,
unsigned int height, unsigned int width) {
using T = Matrix2x2<std::complex<float>>;
- Array4D<T> aterms =
- proxy.allocate_array4d<T>(nr_timeslots, nr_stations, height, width);
-
- for (unsigned t = 0; t < nr_timeslots; t++) {
- for (unsigned ant = 0; ant < nr_stations; ant++) {
- for (unsigned y = 0; y < height; y++) {
- for (unsigned x = 0; x < width; x++) {
- float scale = ((float)(t + 1) / nr_timeslots);
- std::complex<float> valueXX = std::complex<float>(scale * 1.0, 1.1);
- std::complex<float> valueXY = std::complex<float>(scale * 0.8, 0.9);
- std::complex<float> valueYX = std::complex<float>(scale * 0.6, 1.7);
- std::complex<float> valueYY = std::complex<float>(scale * 0.4, 0.5);
- const Matrix2x2<std::complex<float>> aterm = {valueXX, valueXY,
- valueYX, valueYY};
- aterms(t, ant, y, x) = aterm;
- }
- }
- }
- }
-
+ auto aterms =
+ proxy.allocate_span<T, 4>({nr_timeslots, nr_stations, height, width});
+ init_example_aterms(aterms);
return aterms;
}
-Array1D<unsigned int> get_example_aterm_offsets(proxy::Proxy& proxy,
- unsigned int nr_timeslots,
- unsigned int nr_timesteps) {
+aocommon::xt::Span<unsigned int, 1> get_example_aterm_offsets(
+ proxy::Proxy& proxy, unsigned int nr_timeslots, unsigned int nr_timesteps) {
using T = unsigned int;
- Array1D<T> aterm_offsets = proxy.allocate_array1d<T>(nr_timeslots + 1);
-
- for (unsigned time = 0; time < nr_timeslots; time++) {
- aterm_offsets(time) = time * (nr_timesteps / nr_timeslots);
- }
-
- aterm_offsets(nr_timeslots) = nr_timesteps;
-
+ auto aterm_offsets = proxy.allocate_span<T, 1>({nr_timeslots + 1});
+ init_example_aterm_offsets(aterm_offsets, nr_timesteps);
return aterm_offsets;
}
-Array2D<float> get_example_spheroidal(proxy::Proxy& proxy, unsigned int height,
- unsigned int width) {
+aocommon::xt::Span<float, 2> get_identity_spheroidal(proxy::Proxy& proxy,
+ unsigned int height,
+ unsigned int width) {
using T = float;
- Array2D<T> spheroidal = proxy.allocate_array2d<T>(height, width);
-
- // Evaluate rows
- float y[height];
- for (unsigned i = 0; i < height; i++) {
- float tmp = fabs(-1 + i * 2.0f / float(height));
- y[i] = evaluate_spheroidal(tmp);
- }
-
- // Evaluate columns
- float x[width];
- for (unsigned i = 0; i < width; i++) {
- float tmp = fabs(-1 + i * 2.0f / float(width));
- x[i] = evaluate_spheroidal(tmp);
- }
-
- // Set values
- for (unsigned i = 0; i < height; i++) {
- for (unsigned j = 0; j < width; j++) {
- spheroidal(i, j) = y[i] * x[j];
- }
- }
+ auto spheroidal = proxy.allocate_span<T, 2>({height, width});
+ init_identity_spheroidal(spheroidal);
+ return spheroidal;
+}
+aocommon::xt::Span<float, 2> get_example_spheroidal(proxy::Proxy& proxy,
+ unsigned int height,
+ unsigned int width) {
+ using T = float;
+ auto spheroidal = proxy.allocate_span<T, 2>({height, width});
+ init_example_spheroidal(spheroidal);
return spheroidal;
}
/*
* Default memory allocation
*/
-Array1D<float> get_example_frequencies(unsigned int nr_channels,
- float start_frequency,
- float frequency_increment) {
- Array1D<float> frequencies(nr_channels);
+void init_example_frequencies(aocommon::xt::Span<float, 1>& frequencies,
+ float start_frequency,
+ float frequency_increment) {
+ const size_t nr_channels = frequencies.size();
for (unsigned chan = 0; chan < nr_channels; chan++) {
frequencies(chan) = start_frequency + frequency_increment * chan;
}
-
- return frequencies;
}
-Array4D<std::complex<float>> get_dummy_visibilities(
- unsigned int nr_baselines, unsigned int nr_timesteps,
- unsigned int nr_channels, unsigned int nr_correlations) {
+void init_example_visibilities(
+ aocommon::xt::Span<std::complex<float>, 4>& visibilities,
+ aocommon::xt::Span<UVW<float>, 2>& uvw,
+ aocommon::xt::Span<float, 1>& frequencies, float image_size,
+ unsigned int grid_size, unsigned int nr_point_sources,
+ unsigned int max_pixel_offset, unsigned int random_seed, float amplitude) {
+ const size_t nr_baselines = visibilities.shape(0);
+ const size_t nr_timesteps = visibilities.shape(1);
+ const size_t nr_channels = visibilities.shape(2);
+ const size_t nr_correlations = visibilities.shape(3);
assert(nr_correlations == 2 || nr_correlations == 4);
+ assert(uvw.shape(0) == nr_baselines);
+ assert(uvw.shape(1) == nr_timesteps);
+ assert(frequencies.size() == nr_channels);
- Array4D<std::complex<float>> visibilities(nr_baselines, nr_timesteps,
- nr_channels, nr_correlations);
-
-// Set all visibilities
-#pragma omp parallel for
- for (unsigned bl = 0; bl < nr_baselines; bl++) {
- for (unsigned time = 0; time < nr_timesteps; time++) {
- for (unsigned chan = 0; chan < nr_channels; chan++) {
- if (nr_correlations == 2) {
- visibilities(bl, time, chan, 0) = 1.0f;
- visibilities(bl, time, chan, 1) = 1.0f;
- } else if (nr_correlations == 4) {
- visibilities(bl, time, chan, 0) = 1.0f;
- visibilities(bl, time, chan, 1) = 0.0f;
- visibilities(bl, time, chan, 2) = 0.0f;
- visibilities(bl, time, chan, 3) = 1.0f;
- }
- }
- }
+ if (max_pixel_offset == 0) {
+ max_pixel_offset = grid_size / 2;
}
- return visibilities;
-}
-
-Array4D<std::complex<float>> get_example_visibilities(
- Array2D<UVW<float>>& uvw, Array1D<float>& frequencies, float image_size,
- unsigned int grid_size, unsigned int nr_correlations,
- unsigned int nr_point_sources, unsigned int max_pixel_offset,
- unsigned int random_seed, float amplitude) {
- unsigned int nr_baselines = uvw.get_y_dim();
- unsigned int nr_timesteps = uvw.get_x_dim();
- unsigned int nr_channels = frequencies.get_x_dim();
-
- Array4D<std::complex<float>> visibilities(nr_baselines, nr_timesteps,
- nr_channels, nr_correlations);
- std::fill_n(visibilities.data(), visibilities.size(),
- std::complex<float>{0.0, 0.0});
-
srand(random_seed);
- for (unsigned i = 0; i < nr_point_sources; i++) {
- float x_offset = (random() * (max_pixel_offset)) - (max_pixel_offset / 2);
- float y_offset = (random() * (max_pixel_offset)) - (max_pixel_offset / 2);
+ for (size_t i = 0; i < nr_point_sources; i++) {
+ float x_offset = static_cast<float>(random()) / RAND_MAX;
+ float y_offset = static_cast<float>(random()) / RAND_MAX;
+ x_offset = (x_offset * (max_pixel_offset)) - (max_pixel_offset / 2);
+ y_offset = (y_offset * (max_pixel_offset)) - (max_pixel_offset / 2);
add_pt_src(visibilities, uvw, frequencies, image_size, grid_size, x_offset,
y_offset, amplitude);
}
+}
- return visibilities;
+void init_dummy_visibilities(
+ aocommon::xt::Span<std::complex<float>, 4>& visibilities) {
+ const size_t nr_baselines = visibilities.shape(0);
+ const size_t nr_correlations = visibilities.shape(3);
+ assert(nr_correlations == 2 || nr_correlations == 4);
+
+// Set all visibilities
+#pragma omp parallel for
+ for (size_t bl = 0; bl < nr_baselines; bl++) {
+ if (nr_correlations == 2) {
+ xt::view(visibilities, bl, xt::all(), xt::all(), xt::all()) = 1.0f;
+ } else if (nr_correlations == 4) {
+ xt::view(visibilities, bl, xt::all(), xt::all(), 0) = 1.0f;
+ xt::view(visibilities, bl, xt::all(), xt::all(), 1) = 0.0f;
+ xt::view(visibilities, bl, xt::all(), xt::all(), 2) = 0.0f;
+ xt::view(visibilities, bl, xt::all(), xt::all(), 3) = 1.0f;
+ }
+ }
}
-Array1D<std::pair<unsigned int, unsigned int>> get_example_baselines(
+xt::xtensor<std::pair<unsigned int, unsigned int>, 1> get_example_baselines(
unsigned int nr_stations, unsigned int nr_baselines) {
- Array1D<std::pair<unsigned int, unsigned int>> baselines(nr_baselines);
+ const std::array<size_t, 1> shape{nr_baselines};
+ xt::xtensor<std::pair<unsigned int, unsigned int>, 1> baselines(shape);
- unsigned bl = 0;
+ size_t bl = 0;
- for (unsigned station1 = 0; station1 < nr_stations; station1++) {
- for (unsigned station2 = station1 + 1; station2 < nr_stations; station2++) {
+ for (size_t station1 = 0; station1 < nr_stations; station1++) {
+ for (size_t station2 = station1 + 1; station2 < nr_stations; station2++) {
if (bl >= nr_baselines) {
break;
}
@@ -400,7 +235,7 @@ Data get_example_data(unsigned int max_nr_baselines, unsigned int grid_size,
Data data(layout_file);
// Determine the max baseline length for given grid_size
- float max_uv = data.compute_max_uv(grid_size, nr_channels);
+ const float max_uv = data.compute_max_uv(grid_size, nr_channels);
// Select only baselines up to max_uv meters long
data.limit_max_baseline_length(max_uv);
@@ -412,122 +247,115 @@ Data get_example_data(unsigned int max_nr_baselines, unsigned int grid_size,
return data;
}
-Array4D<Matrix2x2<std::complex<float>>> get_identity_aterms(
- unsigned int nr_timeslots, unsigned int nr_stations, unsigned int height,
- unsigned int width) {
- Array4D<Matrix2x2<std::complex<float>>> aterms(nr_timeslots, nr_stations,
- height, width);
- const Matrix2x2<std::complex<float>> aterm = {1.0f, 0.0f, 0.0f, 1.0f};
-
- for (unsigned t = 0; t < nr_timeslots; t++) {
- for (unsigned ant = 0; ant < nr_stations; ant++) {
- for (unsigned y = 0; y < height; y++) {
- for (unsigned x = 0; x < width; x++) {
- aterms(t, ant, y, x) = aterm;
- }
+void add_pt_src(aocommon::xt::Span<std::complex<float>, 4>& visibilities,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<float, 1>& frequencies,
+ float image_size, unsigned int grid_size, float offset_x,
+ float offset_y, float amplitude) {
+ const size_t nr_baselines = visibilities.shape(0);
+ const size_t nr_timesteps = visibilities.shape(1);
+ const size_t nr_channels = visibilities.shape(2);
+
+ const float l = offset_x * image_size / grid_size;
+ const float m = offset_y * image_size / grid_size;
+
+ for (size_t bl = 0; bl < nr_baselines; bl++) {
+ for (size_t t = 0; t < nr_timesteps; t++) {
+ for (size_t c = 0; c < nr_channels; c++) {
+ const double speed_of_light = 299792458.0;
+ const float u = (frequencies(c) / speed_of_light) * uvw(bl, t).u;
+ const float v = (frequencies(c) / speed_of_light) * uvw(bl, t).v;
+ const std::complex<float> value =
+ amplitude *
+ std::exp(std::complex<float>(0, -2 * M_PI * (u * l + v * m)));
+ xt::view(visibilities, bl, t, c, xt::all()) += value;
}
}
}
+}
- return aterms;
+void init_identity_aterms(
+ aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4>& aterms) {
+ const std::complex<float> valueXX{1.0f, 0.0f};
+ const std::complex<float> valueXY{0.0f, 0.0f};
+ const std::complex<float> valueYX{0.0f, 0.0f};
+ const std::complex<float> valueYY{1.0f, 0.0f};
+ const Matrix2x2<std::complex<float>> aterm{valueXX, valueXY, valueYX,
+ valueYY};
+ aterms.fill(aterm);
}
-Array4D<Matrix2x2<std::complex<float>>> get_example_aterms(
- unsigned int nr_timeslots, unsigned int nr_stations, unsigned int height,
- unsigned int width) {
- assert(height == width);
- Array4D<Matrix2x2<std::complex<float>>> aterms(nr_timeslots, nr_stations,
- height, width);
- init_example_aterms(aterms);
- return aterms;
+void init_example_aterms(
+ aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4>& aterms) {
+ const size_t nr_timeslots = aterms.shape(0);
+
+ for (size_t t = 0; t < nr_timeslots; t++) {
+ const float scale = ((float)(t + 1) / nr_timeslots);
+ std::complex<float> valueXX{scale * 1.0f, 1.1f};
+ std::complex<float> valueXY{scale * 0.8f, 0.9f};
+ std::complex<float> valueYX{scale * 0.6f, 1.7f};
+ std::complex<float> valueYY{scale * 0.4f, 0.5f};
+ const Matrix2x2<std::complex<float>> aterm{valueXX, valueXY, valueYX,
+ valueYY};
+ xt::view(aterms, t, xt::all(), xt::all(), xt::all()) = aterm;
+ }
}
-Array1D<unsigned int> get_example_aterm_offsets(unsigned int nr_timeslots,
- unsigned int nr_timesteps) {
- Array1D<unsigned int> aterm_offsets(nr_timeslots + 1);
+void init_example_aterm_offsets(
+ aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
+ unsigned int nr_timesteps) {
+ const size_t nr_timeslots = aterm_offsets.shape(0) - 1;
- for (unsigned time = 0; time < nr_timeslots; time++) {
+ for (size_t time = 0; time < nr_timeslots; time++) {
aterm_offsets(time) = time * (nr_timesteps / nr_timeslots);
}
aterm_offsets(nr_timeslots) = nr_timesteps;
-
- return aterm_offsets;
}
-Array2D<float> get_identity_spheroidal(unsigned int height,
- unsigned int width) {
- Array2D<float> spheroidal(height, width);
-
- float value = 1.0;
-
- for (unsigned y = 0; y < height; y++) {
- for (unsigned x = 0; x < width; x++) {
- spheroidal(y, x) = value;
- }
- }
-
- return spheroidal;
+void init_identity_spheroidal(aocommon::xt::Span<float, 2>& spheroidal) {
+ spheroidal.fill(1.0f);
}
-Array2D<float> get_example_spheroidal(unsigned int height, unsigned int width) {
- Array2D<float> spheroidal(height, width);
+void init_example_spheroidal(aocommon::xt::Span<float, 2>& spheroidal) {
+ const size_t height = spheroidal.shape(0);
+ const size_t width = spheroidal.shape(1);
// Evaluate rows
float y[height];
for (unsigned i = 0; i < height; i++) {
float tmp = fabs(-1 + i * 2.0f / float(height));
- y[i] = evaluate_spheroidal(tmp);
+ y[i] = ::evaluate_spheroidal(tmp);
}
// Evaluate columns
float x[width];
for (unsigned i = 0; i < width; i++) {
float tmp = fabs(-1 + i * 2.0f / float(width));
- x[i] = evaluate_spheroidal(tmp);
+ x[i] = ::evaluate_spheroidal(tmp);
}
// Set values
- for (unsigned i = 0; i < height; i++) {
- for (unsigned j = 0; j < width; j++) {
+ for (size_t i = 0; i < height; i++) {
+ for (size_t j = 0; j < width; j++) {
spheroidal(i, j) = y[i] * x[j];
}
}
-
- return spheroidal;
}
-Array1D<float> get_zero_shift() {
- Array1D<float> shift(2);
- shift.zero();
- return shift;
-}
+void init_example_baselines(
+ aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>& baselines,
+ unsigned int nr_stations) {
+ const size_t nr_baselines = baselines.size();
+ size_t bl = 0;
-void add_pt_src(Array4D<std::complex<float>>& visibilities,
- Array2D<UVW<float>>& uvw, Array1D<float>& frequencies,
- float image_size, unsigned int grid_size, float offset_x,
- float offset_y, float amplitude) {
- auto nr_baselines = visibilities.get_w_dim();
- auto nr_timesteps = visibilities.get_z_dim();
- auto nr_channels = visibilities.get_y_dim();
- auto nr_correlations = visibilities.get_x_dim();
-
- float l = offset_x * image_size / grid_size;
- float m = offset_y * image_size / grid_size;
-
- for (unsigned bl = 0; bl < nr_baselines; bl++) {
- for (unsigned t = 0; t < nr_timesteps; t++) {
- for (unsigned c = 0; c < nr_channels; c++) {
- const double speed_of_light = 299792458.0;
- float u = (frequencies(c) / speed_of_light) * uvw(bl, t).u;
- float v = (frequencies(c) / speed_of_light) * uvw(bl, t).v;
- std::complex<float> value =
- amplitude *
- exp(std::complex<float>(0, -2 * M_PI * (u * l + v * m)));
- for (unsigned cor = 0; cor < nr_correlations; cor++) {
- visibilities(bl, t, c, cor) += value;
- }
+ for (size_t station1 = 0; station1 < nr_stations; station1++) {
+ for (size_t station2 = station1 + 1; station2 < nr_stations; station2++) {
+ if (bl >= nr_baselines) {
+ break;
}
+ baselines(bl) = std::pair<unsigned int, unsigned int>(station1, station2);
+ bl++;
}
}
}
diff --git a/idg-util/src/initialize/init.h b/idg-util/src/initialize/init.h
index c92edba57251ca81f53f1c913627e0880c86f1cf..936b89417121ddbbc1cade64573c49970eaed60f 100644
--- a/idg-util/src/initialize/init.h
+++ b/idg-util/src/initialize/init.h
@@ -21,83 +21,103 @@ Data get_example_data(unsigned int max_nr_baselines, unsigned int grid_size,
/*
* Memory-allocation is handled by Proxy
*/
-Array1D<float> get_example_frequencies(
+aocommon::xt::Span<float, 1> get_example_frequencies(
proxy::Proxy& proxy, unsigned int nr_channels,
float start_frequency = Data::start_frequency,
float frequency_increment = Data::frequency_increment);
-Array4D<std::complex<float>> get_dummy_visibilities(
+aocommon::xt::Span<std::complex<float>, 4> get_dummy_visibilities(
proxy::Proxy& proxy, unsigned int nr_baselines, unsigned int nr_timesteps,
unsigned int nr_channels, unsigned int nr_correlations);
-Array4D<std::complex<float>> get_example_visibilities(
- proxy::Proxy& proxy, Array2D<UVW<float>>& uvw, Array1D<float>& frequencies,
- float image_size, unsigned int nr_correlations, unsigned int grid_size,
- unsigned int nr_polarizations, unsigned int nr_point_sources = 4,
- unsigned int max_pixel_offset = -1, unsigned int random_seed = 2,
- float amplitude = 1);
+aocommon::xt::Span<std::complex<float>, 4> get_example_visibilities(
+ proxy::Proxy& proxy, aocommon::xt::Span<UVW<float>, 2>& uvw,
+ aocommon::xt::Span<float, 1>& frequencies, float image_size,
+ unsigned int nr_correlations, unsigned int grid_size,
+ unsigned int nr_point_sources = 4, unsigned int max_pixel_offset = 0,
+ unsigned int random_seed = 2, float amplitude = 1);
-Array1D<std::pair<unsigned int, unsigned int>> get_example_baselines(
- proxy::Proxy& proxy, unsigned int nr_stations, unsigned int nr_baselines);
+aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>
+get_example_baselines(proxy::Proxy& proxy, unsigned int nr_stations,
+ unsigned int nr_baselines);
-Array4D<Matrix2x2<std::complex<float>>> get_identity_aterms(
+aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4> get_identity_aterms(
proxy::Proxy& proxy, unsigned int nr_timeslots, unsigned int nr_stations,
unsigned int height, unsigned int width);
-Array4D<Matrix2x2<std::complex<float>>> get_example_aterms(
+aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4> get_example_aterms(
proxy::Proxy& proxy, unsigned int nr_timeslots, unsigned int nr_stations,
unsigned int height, unsigned int width);
-Array1D<unsigned int> get_example_aterm_offsets(proxy::Proxy& proxy,
- unsigned int nr_timeslots,
- unsigned int nr_timesteps);
+aocommon::xt::Span<unsigned int, 1> get_example_aterm_offsets(
+ proxy::Proxy& proxy, unsigned int nr_timeslots, unsigned int nr_timesteps);
-Array2D<float> get_example_spheroidal(proxy::Proxy& proxy, unsigned int height,
- unsigned int width);
+aocommon::xt::Span<float, 2> get_identity_spheroidal(proxy::Proxy& proxy,
+ unsigned int height,
+ unsigned int width);
-Array1D<float> get_zero_shift();
+aocommon::xt::Span<float, 2> get_example_spheroidal(proxy::Proxy& proxy,
+ unsigned int height,
+ unsigned int width);
/*
* Default memory allocation
*/
-Array1D<float> get_example_frequencies(
- unsigned int nr_channels, float start_frequency = Data::start_frequency,
+void init_example_frequencies(
+ aocommon::xt::Span<float, 1>& frequencies,
+ float start_frequency = Data::start_frequency,
float frequency_increment = Data::frequency_increment);
-Array4D<std::complex<float>> get_dummy_visibilities(
- unsigned int nr_baselines, unsigned int nr_timesteps,
- unsigned int nr_channels, unsigned int nr_correlations);
+void init_example_visibilities(
+ aocommon::xt::Span<std::complex<float>, 4>& visibilities,
+ aocommon::xt::Span<UVW<float>, 2>& uvw,
+ aocommon::xt::Span<float, 1>& frequencies, float image_size,
+ unsigned int grid_size, unsigned int nr_point_sources = 4,
+ unsigned int max_pixel_offset = 0, unsigned int random_seed = 2,
+ float amplitude = 1);
-Array4D<std::complex<float>> get_example_visibilities(
- Array2D<UVW<float>>& uvw, Array1D<float>& frequencies, float image_size,
- unsigned int grid_size, unsigned int nr_correlations,
- unsigned int nr_point_sources = 4, unsigned int max_pixel_offset = -1,
- unsigned int random_seed = 2, float amplitude = 1);
+void init_dummy_visibilities(
+ aocommon::xt::Span<std::complex<float>, 4>& visibilities);
-Array1D<std::pair<unsigned int, unsigned int>> get_example_baselines(
+xt::xtensor<std::pair<unsigned int, unsigned int>, 1> get_example_baselines(
unsigned int nr_stations, unsigned int nr_baselines);
-Array4D<Matrix2x2<std::complex<float>>> get_identity_aterms(
+xt::xtensor<Matrix2x2<std::complex<float>>, 4> get_identity_aterms(
unsigned int nr_timeslots, unsigned int nr_stations, unsigned int height,
unsigned int width);
-Array4D<Matrix2x2<std::complex<float>>> get_example_aterms(
+xt::xtensor<Matrix2x2<std::complex<float>>, 4> get_example_aterms(
unsigned int nr_timeslots, unsigned int nr_stations, unsigned int height,
unsigned int width);
-Array1D<unsigned int> get_example_aterm_offsets(unsigned int nr_timeslots,
- unsigned int nr_timesteps);
+xt::xtensor<unsigned int, 1> get_example_aterm_offsets(
+ unsigned int nr_timeslots, unsigned int nr_timesteps);
-Array2D<float> get_identity_spheroidal(unsigned int height, unsigned int width);
+void init_identity_spheroidal(aocommon::xt::Span<float, 2>& spheroidal);
-Array2D<float> get_example_spheroidal(unsigned int height, unsigned int width);
+void init_example_spheroidal(aocommon::xt::Span<float, 2>& spheroidal);
float evaluate_spheroidal(float nu);
-void add_pt_src(Array4D<std::complex<float>>& visibilities,
- Array2D<UVW<float>>& uvw, Array1D<float>& frequencies,
- float image_size, unsigned int grid_size, float x, float y,
- float amplitude);
+void add_pt_src(aocommon::xt::Span<std::complex<float>, 4>& visibilities,
+ const aocommon::xt::Span<UVW<float>, 2>& uvw,
+ const aocommon::xt::Span<float, 1>& frequencies,
+ float image_size, unsigned int grid_size, float offset_x,
+ float offset_y, float amplitude);
+
+void init_identity_aterms(
+ aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4>& aterms);
+
+void init_example_aterms(
+ aocommon::xt::Span<Matrix2x2<std::complex<float>>, 4>& aterms);
+
+void init_example_aterm_offsets(
+ aocommon::xt::Span<unsigned int, 1>& aterm_offsets,
+ unsigned int nr_timesteps);
+
+void init_example_baselines(
+ aocommon::xt::Span<std::pair<unsigned int, unsigned int>, 1>& baselines,
+ unsigned int nr_stations);
} // namespace idg
diff --git a/idg-util/src/initialize/initc.h b/idg-util/src/initialize/initc.h
index c67119191e64eeabda644abe12bbeee3c489a0e8..3abca012bca9f122795ab8b55e384ba6380b3c50 100644
--- a/idg-util/src/initialize/initc.h
+++ b/idg-util/src/initialize/initc.h
@@ -4,72 +4,100 @@
extern "C" {
void utils_init_identity_spheroidal(void* ptr, int subgrid_size) {
- idg::Array2D<float> spheroidal =
- idg::get_identity_spheroidal(subgrid_size, subgrid_size);
- memcpy(ptr, spheroidal.data(), spheroidal.bytes());
+ const std::array<size_t, 2> spheroidal_shape{
+ static_cast<size_t>(subgrid_size), static_cast<size_t>(subgrid_size)};
+ auto spheroidal =
+ aocommon::xt::CreateSpan(reinterpret_cast<float*>(ptr), spheroidal_shape);
+ idg::init_identity_spheroidal(spheroidal);
}
void utils_init_example_spheroidal(void* ptr, int subgrid_size) {
- idg::Array2D<float> spheroidal =
- idg::get_example_spheroidal(subgrid_size, subgrid_size);
- memcpy(ptr, spheroidal.data(), spheroidal.bytes());
+ const std::array<size_t, 2> spheroidal_shape{
+ static_cast<size_t>(subgrid_size), static_cast<size_t>(subgrid_size)};
+ auto spheroidal =
+ aocommon::xt::CreateSpan(reinterpret_cast<float*>(ptr), spheroidal_shape);
+ idg::init_example_spheroidal(spheroidal);
}
void utils_init_example_frequencies(void* ptr, int nr_channels) {
- idg::Array1D<float> frequencies = idg::get_example_frequencies(nr_channels);
- memcpy(ptr, frequencies.data(), frequencies.bytes());
+ const std::array<size_t, 1> frequencies_shape{
+ static_cast<size_t>(nr_channels)};
+ auto frequencies = aocommon::xt::CreateSpan(reinterpret_cast<float*>(ptr),
+ frequencies_shape);
+ idg::init_example_frequencies(frequencies);
}
void utils_init_dummy_visibilities(void* ptr, int nr_baselines,
int nr_timesteps, int nr_channels,
int nr_correlations) {
- idg::Array4D<std::complex<float>> visibilities = idg::get_dummy_visibilities(
- nr_baselines, nr_timesteps, nr_channels, nr_correlations);
- memcpy(ptr, visibilities.data(), visibilities.bytes());
+ const std::array<size_t, 4> visibilities_shape{
+ static_cast<size_t>(nr_baselines), static_cast<size_t>(nr_timesteps),
+ static_cast<size_t>(nr_channels), static_cast<size_t>(nr_correlations)};
+ auto visibilities = aocommon::xt::CreateSpan(
+ reinterpret_cast<std::complex<float>*>(ptr), visibilities_shape);
+ idg::init_dummy_visibilities(visibilities);
}
void utils_add_pt_src(float x, float y, float amplitude, int nr_baselines,
int nr_timesteps, int nr_channels, int nr_correlations,
float image_size, int grid_size, void* uvw,
void* frequencies, void* visibilities) {
- typedef idg::UVW<float> UVWType;
- idg::Array4D<std::complex<float>> visibilities_(
- reinterpret_cast<std::complex<float>*>(visibilities), nr_baselines,
- nr_timesteps, nr_channels, nr_correlations);
- idg::Array2D<UVWType> uvw_((UVWType*)uvw, nr_baselines, nr_timesteps);
- idg::Array1D<float> frequencies_((float*)frequencies, nr_channels);
- idg::add_pt_src(visibilities_, uvw_, frequencies_, image_size, grid_size, x,
- y, amplitude);
+ const std::array<size_t, 1> frequencies_shape{
+ static_cast<size_t>(nr_channels)};
+ auto frequencies_span = aocommon::xt::CreateSpan(
+ reinterpret_cast<float*>(frequencies), frequencies_shape);
+ const std::array<size_t, 2> uvw_shape{static_cast<size_t>(nr_baselines),
+ static_cast<size_t>(nr_timesteps)};
+ auto uvw_span = aocommon::xt::CreateSpan(
+ reinterpret_cast<idg::UVW<float>*>(uvw), uvw_shape);
+ const std::array<size_t, 4> visibilities_shape{
+ static_cast<size_t>(nr_baselines), static_cast<size_t>(nr_timesteps),
+ static_cast<size_t>(nr_channels), static_cast<size_t>(nr_correlations)};
+ auto visibilities_span = aocommon::xt::CreateSpan(
+ reinterpret_cast<std::complex<float>*>(visibilities), visibilities_shape);
+ idg::add_pt_src(visibilities_span, uvw_span, frequencies_span, image_size,
+ grid_size, x, y, amplitude);
}
void utils_init_identity_aterms(void* ptr, int nr_timeslots, int nr_stations,
int subgrid_size, int nr_correlations) {
- idg::Array4D<idg::Matrix2x2<std::complex<float>>> aterms =
- idg::get_identity_aterms(nr_timeslots, nr_stations, subgrid_size,
- subgrid_size);
- memcpy(ptr, aterms.data(), aterms.bytes());
+ const std::array<size_t, 4> aterms_shape{
+ static_cast<size_t>(nr_timeslots), static_cast<size_t>(nr_stations),
+ static_cast<size_t>(subgrid_size), static_cast<size_t>(subgrid_size)};
+ using T = idg::Matrix2x2<std::complex<float>>;
+ auto aterms =
+ aocommon::xt::CreateSpan(reinterpret_cast<T*>(ptr), aterms_shape);
+ idg::init_identity_aterms(aterms);
}
void utils_init_example_aterms(void* ptr, int nr_timeslots, int nr_stations,
int subgrid_size, int nr_correlations) {
- idg::Array4D<idg::Matrix2x2<std::complex<float>>> aterms =
- idg::get_example_aterms(nr_timeslots, nr_stations, subgrid_size,
- subgrid_size);
- memcpy(ptr, aterms.data(), aterms.bytes());
+ const std::array<size_t, 4> aterms_shape{
+ static_cast<size_t>(nr_timeslots), static_cast<size_t>(nr_stations),
+ static_cast<size_t>(subgrid_size), static_cast<size_t>(subgrid_size)};
+ using Aterm = idg::Matrix2x2<std::complex<float>>;
+ auto aterms =
+ aocommon::xt::CreateSpan(reinterpret_cast<Aterm*>(ptr), aterms_shape);
+ idg::init_example_aterms(aterms);
}
void utils_init_example_aterm_offsets(void* ptr, int nr_timeslots,
int nr_timesteps) {
- idg::Array1D<unsigned int> aterm_offsets =
- idg::get_example_aterm_offsets(nr_timeslots, nr_timesteps);
- memcpy(ptr, aterm_offsets.data(), aterm_offsets.bytes());
+ const std::array<size_t, 1> aterm_offsets_shape{
+ static_cast<size_t>(nr_timeslots) + 1};
+ auto aterm_offsets = aocommon::xt::CreateSpan(
+ reinterpret_cast<unsigned int*>(ptr), aterm_offsets_shape);
+ idg::init_example_aterm_offsets(aterm_offsets, nr_timesteps);
}
void utils_init_example_baselines(void* ptr, int nr_stations,
int nr_baselines) {
- idg::Array1D<std::pair<unsigned int, unsigned int>> baselines =
- idg::get_example_baselines(nr_stations, nr_baselines);
- memcpy(ptr, baselines.data(), baselines.bytes());
+ const std::array<size_t, 1> baselines_shape{
+ static_cast<size_t>(nr_baselines)};
+ auto baselines = aocommon::xt::CreateSpan(
+ reinterpret_cast<std::pair<unsigned int, unsigned int>*>(ptr),
+ baselines_shape);
+ idg::init_example_baselines(baselines, nr_stations);
}
} // end extern "C"
diff --git a/idg-util/src/python/util.py b/idg-util/src/python/util.py
index 6daba72abe31e64ac37d5726f97b7a17cd7e840c..05efbd2d2dc4387b42eee0709517c80c9e74c1a7 100644
--- a/idg-util/src/python/util.py
+++ b/idg-util/src/python/util.py
@@ -18,31 +18,6 @@ from idg.idgtypes import *
lib = ctypes.cdll.LoadLibrary('libidg-util.so')
-def resize_spheroidal(spheroidal, size, dtype=np.float32):
- """
- Resize a spheroidal
-
- :param spheroidal: Input spheroidal
- :type spheroidal: np.arrayd(type=float) (two dimensions)
- :param size: New size along one axis
- :type size: int
- :param dtype: new dtype, defaults to np.float32
- :type dtype: np.dtype, optional
- :return: New spheroidal
- :rtype: np.array(dtype=float) (two dimensions)
- """
- if spheroidal.shape[1] != spheroidal.shape[0]:
- raise ValueError("Input spheroidal size should be square")
- subgrid_size = spheroidal.shape[0]
- tmp = spheroidal.astype(np.float32)
- result = np.zeros(shape=(size, size), dtype=np.float32)
- lib.utils_resize_spheroidal(tmp.ctypes.data_as(ctypes.c_void_p),
- ctypes.c_int(subgrid_size),
- result.ctypes.data_as(ctypes.c_void_p),
- ctypes.c_int(size))
- return result.astype(dtype)
-
-
def nr_baselines_to_nr_stations(nr_baselines):
"""
Convert number of baselines to number of stations, assuming that all station