1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
#include "purify/config.h"
#include "purify/types.h"
#include <ctime>
#include "purify/directories.h"
#include "purify/logging.h"
#include "purify/measurement_operator_factory.h"
#include "purify/utilities.h"
using namespace purify;
using namespace purify::notinstalled;
int main(int nargs, char const **args) {
auto const session = sopt::mpi::init(nargs, args);
auto const world = sopt::mpi::Communicator::World();
purify::logging::set_level("debug");
t_int const conj = (nargs > 1) ? std::stod(static_cast<std::string>(args[1])) : 0;
const t_int iters = (nargs > 2) ? std::stod(static_cast<std::string>(args[2])) : 1;
// Gridding example
auto const oversample_ratio = 2;
auto const Ju = 4;
auto const Jw = 100;
const t_real FoV = 25 * 60 * 60;
std::string const results = output_filename("times_" + std::to_string(world.size()) +
((static_cast<bool>(conj)) ? "_conj" : "") + ".txt");
std::ofstream out(world.is_root() ? results : output_filename("empty.txt"));
if (world.is_root()) {
out.precision(13);
}
for (t_int nvis_iters = 6; nvis_iters < 9; nvis_iters++) {
for (t_int w_iter = 1; w_iter < 4; w_iter++) {
for (t_int cell_iters = 8; cell_iters < 13; cell_iters++) {
t_int const number_of_vis = std::pow(10, nvis_iters);
const t_real rms_w = 50. * w_iter; // lambda
t_int const imsizex = std::pow(2, cell_iters);
t_int const imsizey = imsizex;
const t_real cell = FoV / static_cast<t_real>(imsizex);
auto const kernel = "kb";
t_uint const number_of_pixels = imsizex * imsizey;
// Generating random uv(w) coverage
t_real const sigma_m = constant::pi / 3;
auto uv_data = utilities::random_sample_density(std::floor(number_of_vis / world.size()), 0,
sigma_m, rms_w);
if (static_cast<bool>(conj)) uv_data = utilities::conjugate_w(uv_data);
const auto cost = [](t_real x) -> t_real { return std::abs(x * x); };
uv_data = utilities::w_stacking(uv_data, world, 100, cost, 1e-3);
uv_data.units = utilities::vis_units::radians;
const Vector<t_complex> image = Vector<t_complex>::Random(number_of_pixels);
t_real con_time = 0;
t_real app_time = 0;
for (t_int i = 0; i < iters; i++) {
auto begin = std::chrono::high_resolution_clock::now();
const auto measure_op = factory::measurement_operator_factory<Vector<t_complex>>(
factory::distributed_measurement_operator::mpi_distribute_image, uv_data, imsizey,
imsizex, cell, cell, oversample_ratio, kernels::kernel_from_string.at(kernel), Ju, Jw,
true, 1e-6, 1e-6, dde_type::wkernel_radial);
auto end = std::chrono::high_resolution_clock::now();
con_time += world.all_reduce<t_real>(
std::chrono::duration_cast<std::chrono::nanoseconds>(end - begin).count() * 1e-9,
MPI_MAX); // total time for construction to run in seconds
world.barrier();
begin = std::chrono::high_resolution_clock::now();
Vector<t_complex> const measurements =
measure_op->adjoint() * (*measure_op * image).eval();
end = std::chrono::high_resolution_clock::now();
app_time += world.all_reduce<t_real>(
std::chrono::duration_cast<std::chrono::nanoseconds>(end - begin).count() * 1e-9,
MPI_MAX); // total time for application to run in seconds
}
con_time /= static_cast<t_real>(iters);
app_time /= static_cast<t_real>(iters);
if (world.is_root()) {
out << number_of_vis << " " << rms_w << " " << imsizex << " " << con_time << " "
<< app_time << std::endl; // in seconds
PURIFY_LOW_LOG("Mean Construction Time: {}, Mean Application Time: {}", con_time,
app_time);
}
}
}
}
out.close();
}
|
