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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
|
#include "purify/config.h"
#include <iostream>
#include "catch.hpp"
#include "purify/directories.h"
#include "purify/types.h"
#include "purify/logging.h"
#include "purify/kernels.h"
#include "purify/directories.h"
#include "purify/pfitsio.h"
#include "purify/wide_field_utilities.h"
#include "purify/wkernel_integration.h"
using namespace purify;
using namespace purify::notinstalled;
using namespace purify;
int main(int nargs, char const **args) {
auto const output_name = (nargs > 2) ? static_cast<std::string>(args[1]) : "kernel";
#define ARGS_MACRO(NAME, ARGN, VALUE, TYPE) \
auto const NAME = \
static_cast<TYPE>((nargs > ARGN) ? std::stod(static_cast<std::string>(args[ARGN])) : VALUE);
ARGS_MACRO(w, 2, 10, t_real)
ARGS_MACRO(absolute_error, 3, 1e-2, t_real)
ARGS_MACRO(imsize, 4, 1024, t_uint)
ARGS_MACRO(cell, 5, 20, t_real)
ARGS_MACRO(radial, 6, false, bool)
#undef ARGS_MACRO
purify::logging::set_level("debug");
t_uint const J = 4;
t_int const Jw = 30;
const t_real oversample_ratio = 2;
PURIFY_LOW_LOG("image size: {}", imsize);
PURIFY_LOW_LOG("cell size: {}", cell);
const t_real du = widefield::pixel_to_lambda(cell, imsize, oversample_ratio);
PURIFY_LOW_LOG("1/du: {}", 1. / du);
auto const normu = [=](const t_real l) -> t_real { return kernels::ft_kaiser_bessel(l, J); };
const t_uint max_evaluations = 1e8;
const t_real relative_error = 0;
t_uint e;
const t_real norm =
(not radial)
? std::abs(projection_kernels::exact_w_projection_integration(
0, 0, 0, du, du, oversample_ratio, normu, normu, 1e9, 0, 1e-8,
integration::method::h, e))
: std::abs(projection_kernels::exact_w_projection_integration_1d(
0, 0, 0, du, oversample_ratio, normu, 1e9, 0, 1e-7, integration::method::h, e));
auto const ftkernelu = [=](const t_real l) -> t_real {
return kernels::ft_kaiser_bessel(l, J) / norm;
};
auto const ftkernelv = [=](const t_real l) -> t_real { return kernels::ft_kaiser_bessel(l, J); };
t_uint total = 0;
t_uint rtotal = 0;
t_int const Ju = 20;
t_real const upsample = 5;
t_int const Ju_max = std::floor(Ju * upsample);
const t_int Jw_max =
100; // std::min<int>(std::max<int>(std::floor(std::abs(w * 1./du)), J), Jw);
Vector<t_complex> kernel = Vector<t_complex>::Zero(Ju_max * Jw_max);
Vector<t_real> evals = Vector<t_real>::Zero(Ju_max * Jw_max);
Vector<t_real> support = Vector<t_real>::Zero(Jw_max);
PURIFY_LOW_LOG("w: {}", w);
PURIFY_LOW_LOG("Kernel size: {} x {}", Ju_max, Jw_max);
PURIFY_LOW_LOG("FoV (deg): {}", imsize * cell / 60 / 60);
PURIFY_LOW_LOG("du: {}", du);
#pragma omp parallel for collapse(2)
for (int i = 0; i < Jw_max; i++) {
for (int j = 0; j < Ju_max; j++) {
t_uint evaluations = 0;
t_uint revaluations = 0;
t_uint e;
kernel(j * Jw_max + i) =
(not radial)
? projection_kernels::exact_w_projection_integration(
j / upsample, 0, w * i / Jw_max, du, du, oversample_ratio, ftkernelu, ftkernelv,
max_evaluations, absolute_error, 0, integration::method::h, evaluations)
: projection_kernels::exact_w_projection_integration_1d(
j / upsample, 0, w * i / Jw_max, du, oversample_ratio, ftkernelu,
max_evaluations, absolute_error, 0, integration::method::p, evaluations);
evals(j * Jw_max + i) = evaluations;
#pragma omp critical(print)
{
PURIFY_LOW_LOG("u : {}, w : {}", j / upsample, w * i / Jw_max);
PURIFY_LOW_LOG("Kernel value: {}", kernel(j * Jw_max + i));
PURIFY_LOW_LOG("Evaluations: {}", evaluations);
}
if (kernel(j * Jw_max + i) != kernel(j * Jw_max + i))
throw std::runtime_error("Kernel value is a nan.");
total += evaluations;
}
}
PURIFY_LOW_LOG("Total Evaluations: {}", total);
PURIFY_LOW_LOG("FoV (deg): {}", imsize * cell / 60 / 60);
PURIFY_LOW_LOG("du: {}", du);
pfitsio::write2d(kernel.real(), Jw_max, Ju_max, output_name + "_real.fits");
pfitsio::write2d(kernel.imag(), Jw_max, Ju_max, output_name + "_imag.fits");
pfitsio::write2d(kernel.cwiseAbs(), Jw_max, Ju_max, output_name + ".fits");
pfitsio::write2d(evals, Jw_max, Ju_max, output_name + "_evals.fits");
}
|
