From da70313ca58118820332bb47de963ce9d1516849 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Sun, 10 Feb 2019 09:13:29 +0000 Subject: [PATCH 01/22] added basic RM operator --- cpp/purify/CMakeLists.txt | 3 +- cpp/purify/rm_operators.cc | 47 ++++++++ cpp/purify/rm_operators.h | 216 +++++++++++++++++++++++++++++++++++++ 3 files changed, 265 insertions(+), 1 deletion(-) create mode 100644 cpp/purify/rm_operators.cc create mode 100644 cpp/purify/rm_operators.h diff --git a/cpp/purify/CMakeLists.txt b/cpp/purify/CMakeLists.txt index a21a95efb..03a936a1c 100644 --- a/cpp/purify/CMakeLists.txt +++ b/cpp/purify/CMakeLists.txt @@ -28,12 +28,13 @@ set(HEADERS wide_field_utilities.h wkernel_integration.h wproj_operators.h + rm_operators.h "${PROJECT_BINARY_DIR}/include/purify/config.h") set(SOURCES utilities.cc pfitsio.cc kernels.cc wproj_utilities.cc operators.cc uvfits.cc yaml-parser.cc read_measurements.cc distribute.cc integration.cc wide_field_utilities.cc wkernel_integration.cc - wproj_operators.cc) + wproj_operators.cc rm_operators.cc) if(TARGET casacore::ms) list(APPEND SOURCES casacore.cc) diff --git a/cpp/purify/rm_operators.cc b/cpp/purify/rm_operators.cc new file mode 100644 index 000000000..81665b1ba --- /dev/null +++ b/cpp/purify/rm_operators.cc @@ -0,0 +1,47 @@ +#include "purify/rm_operators.h" +namespace purify { + +namespace details { +//! Construct gridding matrix +Sparse init_gridding_matrix_1d(const Vector &u, const Vector &weights, + const t_uint imsizex_, const t_real oversample_ratio, + const std::function kernelu, + const t_uint Ju) { + const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio); + const t_uint rows = u.size(); + const t_uint cols = ftsizeu_; + + Sparse interpolation_matrix(rows, cols); + interpolation_matrix.reserve(Vector::Constant(rows, Ju)); + + const t_complex I(0, 1); + const t_int ju_max = std::min(Ju, ftsizeu_); + // If I collapse this for loop there is a crash when using MPI... Sparse<>::insert() doesn't work + // right +#pragma omp parallel for + for (t_int m = 0; m < rows; ++m) { + for (t_int ju = 1; ju < ju_max + 1; ++ju) { + const t_real k_u = std::floor(u(m) - ju_max * 0.5); + const t_uint q = utilities::mod(k_u + ju, ftsizeu_); + const t_uint index = utilities::sub2ind(0, q, 1, ftsizeu_); + interpolation_matrix.insert(m, index) = std::exp(-2 * constant::pi * I * (k_u + ju) * 0.5) * + kernelu(u(m) - (k_u + ju)) * weights(m); + } + } + return interpolation_matrix; +} + +Image init_correction1d(const t_real oversample_ratio, const t_uint imsizex_, + const std::function ftkernelu) { + const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio); + const t_uint x_start = std::floor(ftsizeu_ * 0.5 - imsizex_ * 0.5); + + Array range; + range.setLinSpaced(ftsizeu_, 0.5, ftsizeu_ - 0.5); + return (1e0 / range.segment(x_start, imsizex_).unaryExpr([&](t_real x) { + return ftkernelu(x / static_cast(ftsizeu_) - 0.5); + })) * + t_complex(1., 0.); +} +} // namespace details +} // namespace purify diff --git a/cpp/purify/rm_operators.h b/cpp/purify/rm_operators.h new file mode 100644 index 000000000..d3aaf4aca --- /dev/null +++ b/cpp/purify/rm_operators.h @@ -0,0 +1,216 @@ +#ifndef PURIFY_RM_OPERATOR_H +#define PURIFY_RM_OPERATOR_H + +#include "purify/config.h" +#include "purify/types.h" + +#include "purify/operators.h" + +namespace purify { + +namespace details { + +//! Construct gridding matrix +Sparse init_gridding_matrix_1d(const Vector &u, const Vector &weights, + const t_uint imsizex_, const t_real oversample_ratio, + const std::function kernelu, + const t_uint Ju); + +//! Given the Fourier transform of a gridding kernel, creates the scaling image for gridding +//! correction. +Image init_correction1d(const t_real oversample_ratio, const t_uint imsizex_, + const std::function ftkernelu); + +} // namespace details + +namespace operators { + +//! constructs lambdas that apply degridding matrix with adjoint +template +std::tuple, sopt::OperatorFunction> init_gridding_matrix_1d( + ARGS &&... args) { + const std::shared_ptr> interpolation_matrix = + std::make_shared>( + details::init_gridding_matrix_1d(std::forward(args)...)); + const std::shared_ptr> adjoint = + std::make_shared>(interpolation_matrix->adjoint()); + + return std::make_tuple( + [=](T &output, const T &input) { + output = utilities::sparse_multiply_matrix(*interpolation_matrix, input); + }, + [=](T &output, const T &input) { + output = utilities::sparse_multiply_matrix(*adjoint, input); + }); +} + +//! Construsts zero padding operator +template +std::tuple, sopt::OperatorFunction> init_zero_padding_1d( + const Image &S, const t_real oversample_ratio) { + const t_uint imsizex_ = S.cols(); + const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio); + const t_uint x_start = std::floor(ftsizeu_ * 0.5 - imsizex_ * 0.5); + auto direct = [=](T &output, const T &x) { + assert(x.size() == imsizex_); + output = Vector::Zero(ftsizeu_); +#pragma omp parallel for + for (t_uint j = 0; j < imsizex_; j++) output(x_start + j) = S(j) * x(j); + }; + auto indirect = [=](T &output, const T &x) { + assert(x.size() == ftsizeu_); + output = T::Zero(imsizex_); +#pragma omp parallel for + for (t_uint j = 0; j < imsizex_; j++) output(j) = std::conj(S(j)) * x(x_start + j); + }; + return std::make_tuple(direct, indirect); +} // namespace operators + +//! Construsts FFT operator +template +std::tuple, sopt::OperatorFunction> init_FFT_1d( + const t_uint imsizex_, const t_real oversample_factor_, + const fftw_plan fftw_plan_flag_ = fftw_plan::measure) { + t_int const ftsizeu_ = std::floor(imsizex_ * oversample_factor_); + t_int plan_flag = (FFTW_MEASURE | FFTW_PRESERVE_INPUT); + switch (fftw_plan_flag_) { + case (fftw_plan::measure): + plan_flag = (FFTW_MEASURE | FFTW_PRESERVE_INPUT); + break; + case (fftw_plan::estimate): + plan_flag = (FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); + break; + } + +#ifdef PURIFY_OPENMP_FFTW + PURIFY_LOW_LOG("Using OpenMP threading with FFTW."); + fftw_init_threads(); +#endif + Vector src = Vector::Zero(ftsizeu_); + Vector dst = Vector::Zero(ftsizeu_); + // creating plans + const auto del = [](fftw_plan_s *plan) { fftw_destroy_plan(plan); }; + // fftw plan with threads needs to be used before each fftw_plan is created +#ifdef PURIFY_OPENMP_FFTW + fftw_plan_with_nthreads(omp_get_max_threads()); +#endif + const std::shared_ptr m_plan_forward( + fftw_plan_dft_1d(ftsizeu_, reinterpret_cast(src.data()), + reinterpret_cast(dst.data()), FFTW_FORWARD, plan_flag), + del); + // fftw plan with threads needs to be used before each fftw_plan is created +#ifdef PURIFY_OPENMP_FFTW + fftw_plan_with_nthreads(omp_get_max_threads()); +#endif + const std::shared_ptr m_plan_inverse( + fftw_plan_dft_1d(ftsizeu_, reinterpret_cast(src.data()), + reinterpret_cast(dst.data()), FFTW_BACKWARD, plan_flag), + del); + auto const direct = [m_plan_forward, ftsizeu_](T &output, const T &input) { + assert(input.size() == ftsizeu_); + output = Vector::Zero(input.size()); + fftw_execute_dft( + m_plan_forward.get(), + const_cast(reinterpret_cast(input.data())), + reinterpret_cast(output.data())); + output /= std::sqrt(output.size()); + }; + auto const indirect = [m_plan_inverse, ftsizeu_](T &output, const T &input) { + assert(input.size() == ftsizeu_); + output = Vector::Zero(input.size()); + fftw_execute_dft( + m_plan_inverse.get(), + const_cast(reinterpret_cast(input.data())), + reinterpret_cast(output.data())); + output /= std::sqrt(output.size()); + }; + return std::make_tuple(direct, indirect); +} + +template +std::tuple, sopt::OperatorFunction> base_padding_and_FFT_1d( + const std::function ftkernelu, const t_uint imsizex, + const t_real oversample_ratio = 2, const fftw_plan &ft_plan = fftw_plan::measure) { + sopt::OperatorFunction directZ, indirectZ; + sopt::OperatorFunction directFFT, indirectFFT; + + const Image S = + purify::details::init_correction1d(oversample_ratio, imsizex, ftkernelu); + PURIFY_LOW_LOG("Building Measurement Operator: WGFZDB"); + PURIFY_LOW_LOG( + "Constructing Zero Padding " + "and Correction Operator: " + "ZDB"); + PURIFY_MEDIUM_LOG("Faraday Depth size (width): {}", imsizex); + PURIFY_MEDIUM_LOG("Oversampling Factor: {}", oversample_ratio); + std::tie(directZ, indirectZ) = purify::operators::init_zero_padding_1d(S, oversample_ratio); + PURIFY_LOW_LOG("Constructing FFT operator: F"); + switch (ft_plan) { + case fftw_plan::measure: + PURIFY_MEDIUM_LOG("Measuring Plans..."); + break; + case fftw_plan::estimate: + PURIFY_MEDIUM_LOG("Estimating Plans..."); + break; + } + std::tie(directFFT, indirectFFT) = + purify::operators::init_FFT_1d(imsizex, oversample_ratio, ft_plan); + auto direct = sopt::chained_operators(directFFT, directZ); + auto indirect = sopt::chained_operators(indirectZ, indirectFFT); + return std::make_tuple(direct, indirect); +} + +template +std::tuple, sopt::OperatorFunction> base_degrid_operator_1d( + const Vector &u, const Vector &weights, const t_uint imsizey, + const t_uint imsizex, const t_real oversample_ratio = 2, + const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4, const t_uint Jv = 4, + const fftw_plan ft_plan = fftw_plan::measure) { + std::function kernelu, ftkernelu; + std::tie(kernelu, ftkernelu) = purify::create_radial_ftkernel(kernel, Ju, oversample_ratio); + sopt::OperatorFunction directFZ, indirectFZ; + std::tie(directFZ, indirectFZ) = + base_padding_and_FFT_1d(ftkernelu, imsizex, oversample_ratio, ft_plan); + sopt::OperatorFunction directG, indirectG; + PURIFY_LOW_LOG("Constructing Weighting and Gridding Operators: WG"); + PURIFY_MEDIUM_LOG("Number of channels: {}", u.size()); + std::tie(directG, indirectG) = purify::operators::init_gridding_matrix_1d( + u, weights, imsizex, oversample_ratio, kernelu, Ju); + auto direct = sopt::chained_operators(directG, directFZ); + auto indirect = sopt::chained_operators(indirectFZ, indirectG); + PURIFY_LOW_LOG("Finished consturction of RM Φ."); + return std::make_tuple(direct, indirect); +} + +} // namespace operators + +namespace measurementoperator { + +//! Returns linear transform that is the standard degridding operator +template +std::shared_ptr> init_degrid_operator_1d( + const Vector &u, const Vector &weights, const t_uint imsizex, + const t_real oversample_ratio = 2, const kernels::kernel kernel = kernels::kernel::kb, + const t_uint Ju = 4) { + const operators::fftw_plan ft_plan = operators::fftw_plan::measure; + std::array N = {0, 1, static_cast(imsizex)}; + std::array M = {0, 1, static_cast(u.size())}; + sopt::OperatorFunction directDegrid, indirectDegrid; + std::tie(directDegrid, indirectDegrid) = purify::operators::base_degrid_operator_1d( + u, weights, imsizex, oversample_ratio, kernel, Ju); + return std::make_shared>(directDegrid, M, indirectDegrid, N); +} + +template +std::shared_ptr> init_degrid_operator_1d( + const utilities::vis_params &uv_vis_input, const t_uint imsizex, + const t_real oversample_ratio = 2, const kernels::kernel kernel = kernels::kernel::kb, + const t_uint Ju = 4) { + auto uv_vis = uv_vis_input; + return init_degrid_operator_1d(uv_vis.u, uv_vis.weights, imsizex, oversample_ratio, kernel, + Ju); +} + +} // namespace measurementoperator +}; // namespace purify +#endif From d1fcf0e2e0fbb6c5635a9d570db6322ba709fbc9 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Sun, 10 Feb 2019 09:15:43 +0000 Subject: [PATCH 02/22] added comment --- cpp/purify/rm_operators.h | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/purify/rm_operators.h b/cpp/purify/rm_operators.h index d3aaf4aca..815b5bb59 100644 --- a/cpp/purify/rm_operators.h +++ b/cpp/purify/rm_operators.h @@ -207,6 +207,7 @@ std::shared_ptr> init_degrid_operator_1d( const t_real oversample_ratio = 2, const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4) { auto uv_vis = uv_vis_input; + //need to find a way to set RM resolution, i.e. Faraday pixel size return init_degrid_operator_1d(uv_vis.u, uv_vis.weights, imsizex, oversample_ratio, kernel, Ju); } From 3ed245c8dfbe4b5c0e4a51f6ec707d536f243aab Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Tue, 23 Apr 2019 12:16:27 +0100 Subject: [PATCH 03/22] fixed compile errors and bugs so tests pass --- cpp/purify/rm_operators.cc | 3 +-- cpp/purify/rm_operators.h | 15 +++++++-------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/cpp/purify/rm_operators.cc b/cpp/purify/rm_operators.cc index 81665b1ba..6a8d00cc4 100644 --- a/cpp/purify/rm_operators.cc +++ b/cpp/purify/rm_operators.cc @@ -22,8 +22,7 @@ Sparse init_gridding_matrix_1d(const Vector &u, const Vector< for (t_int m = 0; m < rows; ++m) { for (t_int ju = 1; ju < ju_max + 1; ++ju) { const t_real k_u = std::floor(u(m) - ju_max * 0.5); - const t_uint q = utilities::mod(k_u + ju, ftsizeu_); - const t_uint index = utilities::sub2ind(0, q, 1, ftsizeu_); + const t_uint index = utilities::mod(k_u + ju, ftsizeu_); interpolation_matrix.insert(m, index) = std::exp(-2 * constant::pi * I * (k_u + ju) * 0.5) * kernelu(u(m) - (k_u + ju)) * weights(m); } diff --git a/cpp/purify/rm_operators.h b/cpp/purify/rm_operators.h index 815b5bb59..d051ce457 100644 --- a/cpp/purify/rm_operators.h +++ b/cpp/purify/rm_operators.h @@ -48,7 +48,7 @@ std::tuple, sopt::OperatorFunction> init_gridding_m template std::tuple, sopt::OperatorFunction> init_zero_padding_1d( const Image &S, const t_real oversample_ratio) { - const t_uint imsizex_ = S.cols(); + const t_uint imsizex_ = S.size(); const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio); const t_uint x_start = std::floor(ftsizeu_ * 0.5 - imsizex_ * 0.5); auto direct = [=](T &output, const T &x) { @@ -162,12 +162,11 @@ std::tuple, sopt::OperatorFunction> base_padding_an template std::tuple, sopt::OperatorFunction> base_degrid_operator_1d( - const Vector &u, const Vector &weights, const t_uint imsizey, - const t_uint imsizex, const t_real oversample_ratio = 2, - const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4, const t_uint Jv = 4, - const fftw_plan ft_plan = fftw_plan::measure) { + const Vector &u, const Vector &weights, const t_uint imsizex, + const t_real oversample_ratio = 2, const kernels::kernel kernel = kernels::kernel::kb, + const t_uint Ju = 4, const fftw_plan ft_plan = fftw_plan::measure) { std::function kernelu, ftkernelu; - std::tie(kernelu, ftkernelu) = purify::create_radial_ftkernel(kernel, Ju, oversample_ratio); + std::tie(ftkernelu, kernelu) = purify::create_radial_ftkernel(kernel, Ju, oversample_ratio); sopt::OperatorFunction directFZ, indirectFZ; std::tie(directFZ, indirectFZ) = base_padding_and_FFT_1d(ftkernelu, imsizex, oversample_ratio, ft_plan); @@ -197,7 +196,7 @@ std::shared_ptr> init_degrid_operator_1d( std::array M = {0, 1, static_cast(u.size())}; sopt::OperatorFunction directDegrid, indirectDegrid; std::tie(directDegrid, indirectDegrid) = purify::operators::base_degrid_operator_1d( - u, weights, imsizex, oversample_ratio, kernel, Ju); + u, weights, imsizex, oversample_ratio, kernel, Ju, ft_plan); return std::make_shared>(directDegrid, M, indirectDegrid, N); } @@ -207,7 +206,7 @@ std::shared_ptr> init_degrid_operator_1d( const t_real oversample_ratio = 2, const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4) { auto uv_vis = uv_vis_input; - //need to find a way to set RM resolution, i.e. Faraday pixel size + // need to find a way to set RM resolution, i.e. Faraday pixel size return init_degrid_operator_1d(uv_vis.u, uv_vis.weights, imsizex, oversample_ratio, kernel, Ju); } From 23f1fa549f53fc0716e5c5f90fe9b2516474253e Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Tue, 23 Apr 2019 12:17:25 +0100 Subject: [PATCH 04/22] added rm operator tests --- cpp/tests/CMakeLists.txt | 1 + cpp/tests/rm_operator.cc | 260 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 261 insertions(+) create mode 100644 cpp/tests/rm_operator.cc diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 1aa2c424d..b43b0ed84 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -19,6 +19,7 @@ include_directories("${PROJECT_SOURCE_DIR}/cpp" "${CMAKE_CURRENT_BINARY_DIR}/inc file(MAKE_DIRECTORY "${PROJECT_BINARY_DIR}/outputs") add_catch_test(measurement_operator LIBRARIES libpurify) +add_catch_test(rm_operator LIBRARIES libpurify) add_catch_test(purify_fitsio LIBRARIES libpurify) add_catch_test(distribute LIBRARIES libpurify) add_catch_test(utils LIBRARIES libpurify) diff --git a/cpp/tests/rm_operator.cc b/cpp/tests/rm_operator.cc new file mode 100644 index 000000000..17dbb69e0 --- /dev/null +++ b/cpp/tests/rm_operator.cc @@ -0,0 +1,260 @@ +#include +#include "catch.hpp" + +#include "purify/directories.h" +#include "purify/kernels.h" +#include "purify/pfitsio.h" +#include "purify/rm_operators.h" +#include "purify/test_data.h" +#include "purify/utilities.h" +#include + +using namespace purify::notinstalled; +using namespace purify; + +TEST_CASE("1d zeropad") { + const t_int imsize = 8; + const t_real oversample_ratio = 2; + Image S = Vector::Constant(imsize, 2); + auto pad = purify::operators::init_zero_padding_1d>(S, oversample_ratio); + const t_int boundary_size = static_cast(imsize * oversample_ratio * 0.5 - imsize * 0.5); + Vector output; + const Vector input = Vector::Constant(imsize, 5); + std::get<0>(pad)(output, input); + CHECK(output.size() == static_cast(imsize * oversample_ratio)); + CAPTURE(boundary_size); + CAPTURE(input); + CAPTURE(output); + CHECK(output.head(boundary_size).isApprox(Vector::Zero(boundary_size))); + CHECK(output.tail(boundary_size).isApprox(Vector::Zero(boundary_size))); + CHECK(output.segment(boundary_size, imsize).isApprox(Vector::Constant(imsize, 10))); + // apply adjoint + Vector b_output; + std::get<1>(pad)(b_output, output); + CHECK(b_output.size() == static_cast(imsize)); + CHECK(b_output.isApprox(Vector::Constant(imsize, 20))); +} + +TEST_CASE("1d FFT") { + const t_int imsize = 128; + const t_real oversample_factor = 2; + const operators::fftw_plan fftw_plan_flag_ = operators::fftw_plan::measure; + auto fft = operators::init_FFT_1d>(imsize, oversample_factor, fftw_plan_flag_); + SECTION("forward") { + Vector input = Vector::Zero(imsize * oversample_factor); + input(static_cast(0)) = 1.; + Vector output; + std::get<0>(fft)(output, input); + const Vector expected = + Vector::Ones(imsize * oversample_factor) / std::sqrt(imsize * oversample_factor); + CAPTURE(output.head(10)); + CAPTURE(expected.head(10)); + CHECK(output.size() == input.size()); + CHECK(output.isApprox(expected, 1e-6)); + } + SECTION("backward") { + Vector input = Vector::Zero(imsize * oversample_factor); + input(static_cast(0)) = 1.; + Vector output; + std::get<1>(fft)(output, input); + const Vector expected = + Vector::Ones(imsize * oversample_factor) / std::sqrt(imsize * oversample_factor); + CAPTURE(output.head(10)); + CAPTURE(expected.head(10)); + CHECK(output.size() == input.size()); + CHECK(output.isApprox(expected, 1e-6)); + } +} + +TEST_CASE("regression_degrid") { + const t_int imsizex = 256; + Vector u = Vector::Random(10) * imsizex / 2; + Vector input = Vector::Zero(imsizex); + input(static_cast(imsizex * 0.5)) = 1.; + const t_uint M = u.size(); + const t_real oversample_ratio = 2; + const t_int ftsizeu = std::floor(imsizex * oversample_ratio); + const t_uint Ju = 8; + + const Vector y = Vector::Ones(u.size()); + CAPTURE(u); + + SECTION("kb") { + const kernels::kernel kernel = kernels::kernel::kb; + const std::shared_ptr>> measure_op = + measurementoperator::init_degrid_operator_1d>( + u, Vector::Ones(M), imsizex, oversample_ratio, kernel, Ju); + const Vector y_test = *measure_op * input; + CAPTURE(y_test.array() / y.array()); + const t_real max_test = y_test.cwiseAbs().mean(); + CAPTURE(y_test / max_test); + CAPTURE(y); + CHECK((y_test / max_test).isApprox((y), 1e-6)); + } + SECTION("pswf") { + const kernels::kernel kernel = kernels::kernel::pswf; + const auto measure_op = measurementoperator::init_degrid_operator_1d>( + u, Vector::Ones(M), imsizex, oversample_ratio, kernel, 6); + const Vector y_test = *measure_op * input; + CAPTURE(y_test.array() / y.array()); + const t_real max_test = y_test.cwiseAbs().mean(); + CAPTURE(y_test / max_test); + CAPTURE(y); + CHECK((y_test / max_test).isApprox((y), 1e-3)); + } + SECTION("gauss") { + const kernels::kernel kernel = kernels::kernel::gauss; + const auto measure_op = measurementoperator::init_degrid_operator_1d>( + u, Vector::Ones(M), imsizex, oversample_ratio, kernel, Ju); + const Vector y_test = *measure_op * input; + CAPTURE(y_test.array() / y.array()); + const t_real max_test = y_test.cwiseAbs().mean(); + CAPTURE(y_test / max_test); + CAPTURE(y); + CHECK((y_test / max_test).isApprox((y), 1e-4)); + } +} + +TEST_CASE("flux units") { + const t_real oversample_ratio = 2; + Vector u = Vector::Random(10); + const t_uint M = u.size(); + const Vector y = Vector::Ones(u.size()); + SECTION("kb") { + const kernels::kernel kernel = kernels::kernel::kb; + for (auto& J : {4, 5, 6, 7, 8}) { + for (auto& imsize : {128, 256, 512}) { + const auto measure_op = measurementoperator::init_degrid_operator_1d>( + u, Vector::Ones(M), imsize, oversample_ratio, kernel, J); + Vector input = Vector::Zero(imsize); + input(static_cast(imsize * 0.5)) = 1.; + const Vector y_test = + (*measure_op * input).eval() * + std::sqrt(static_cast(input.size()) * oversample_ratio); + CAPTURE(y_test.cwiseAbs().mean()); + CAPTURE(y); + CAPTURE(y_test); + CAPTURE(J); + CAPTURE(imsize) + CHECK(y_test.isApprox(y, 1e-3)); + const Vector psf = + (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M; + CHECK(std::real(psf(static_cast(imsize * 0.5))) == Approx(1.).margin(0.001)); + } + } + } + SECTION("pswf") { + const kernels::kernel kernel = kernels::kernel::pswf; + for (auto& J : {4, 5, 6, 7, 8}) { + for (auto& imsize : {128, 256, 512}) { + if (J != 6) + CHECK_THROWS(measurementoperator::init_degrid_operator_1d>( + u * imsize / 2, Vector::Ones(M), imsize, oversample_ratio, kernel, J)); + else { + const auto measure_op = measurementoperator::init_degrid_operator_1d>( + u, Vector::Ones(M), imsize, oversample_ratio, kernel, J); + Vector input = Vector::Zero(imsize); + input(static_cast(imsize * 0.5)) = 1.; + const Vector y_test = + (*measure_op * input).eval() * + std::sqrt(static_cast(input.size()) * oversample_ratio); + CAPTURE(y_test.cwiseAbs().mean()); + CAPTURE(y); + CAPTURE(y_test); + CAPTURE(J); + CAPTURE(imsize) + CHECK(y_test.isApprox(y, 1e-3)); + const Vector psf = + (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M; + CHECK(std::real(psf(static_cast(imsize * 0.5))) == Approx(1.).margin(0.001)); + } + } + } + } + SECTION("gauss") { + const kernels::kernel kernel = kernels::kernel::gauss; + for (auto& J : {4, 5, 6, 7, 8}) { + for (auto& imsize : {128, 256, 512}) { + const auto measure_op = measurementoperator::init_degrid_operator_1d>( + u, Vector::Ones(M), imsize, oversample_ratio, kernel, J); + Vector input = Vector::Zero(imsize); + input(static_cast(imsize * 0.5)) = 1.; + const Vector y_test = + (*measure_op * input).eval() * + std::sqrt(static_cast(input.size()) * oversample_ratio); + CAPTURE(y_test.cwiseAbs().mean()); + CAPTURE(y); + CAPTURE(y_test); + CAPTURE(J); + CAPTURE(imsize) + CHECK(y_test.isApprox(y, 1e-2)); + const Vector psf = + (measure_op->adjoint() * y) * + std::sqrt(input.size() * oversample_ratio ) / M; + CHECK(std::real(psf(static_cast(imsize * 0.5))) == Approx(1.).margin(0.01)); + } + } + } + SECTION("box") { + const kernels::kernel kernel = kernels::kernel::box; + for (auto& J : {4, 5, 6, 7, 8}) { + for (auto& imsize : {128, 256, 512}) { + const auto measure_op = measurementoperator::init_degrid_operator_1d>( + u, Vector::Ones(M), imsize, oversample_ratio, kernel, J); + Vector input = Vector::Zero(imsize); + input(static_cast(imsize * 0.5)) = 1.; + const Vector y_test = + (*measure_op * input).eval() * + std::sqrt(static_cast(input.size()) * oversample_ratio); + CAPTURE(y_test.cwiseAbs().mean()); + CAPTURE(y); + CAPTURE(y_test); + CAPTURE(J); + CAPTURE(imsize) + CHECK(y_test.isApprox(y, 1e-3)); + const Vector psf = + (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M; + CHECK(std::real(psf(static_cast(imsize * 0.5))) == Approx(1.).margin(0.001)); + } + } + } +} + +TEST_CASE("normed operator") { + const t_real oversample_ratio = 2; + const t_int power_iters = 1000; + const t_real power_tol = 1e-4; + Vector u = Vector::Random(10); + const t_uint M = u.size(); + const Vector y = Vector::Ones(u.size()); + SECTION("kb") { + const kernels::kernel kernel = kernels::kernel::kb; + for (auto& J : {4, 5, 6, 7, 8}) { + for (auto& imsize : {128, 256, 512}) { + const Vector init = Vector::Ones(imsize); + auto power_method_result = sopt::algorithm::normalise_operator>( + measurementoperator::init_degrid_operator_1d>( + u * imsize / 2, Vector::Ones(M), imsize, oversample_ratio, kernel, J), + power_iters, power_tol, init); + const t_real norm = std::get<0>(power_method_result); + const std::shared_ptr>> measure_op = + std::get<2>(power_method_result); + + Vector input = Vector::Zero(imsize); + input(static_cast(imsize * 0.5)) = 1.; + const Vector y_test = + (*measure_op * input).eval() * + std::sqrt(static_cast(input.size()) * oversample_ratio); + CAPTURE(y_test.cwiseAbs().mean()); + CAPTURE(y); + CAPTURE(y_test); + CAPTURE(J); + CAPTURE(imsize) + CHECK((y_test * norm).isApprox(y, 1e-3)); + const Vector psf = + (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M * norm; + CHECK(std::real(psf(static_cast(imsize * 0.5))) == Approx(1.).margin(0.001)); + } + } + } +} From 00ee7741fbdf395a599ca10b7991f4ed98873bc7 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Tue, 23 Apr 2019 18:58:29 +0100 Subject: [PATCH 05/22] added rm integration kernel --- cpp/purify/CMakeLists.txt | 3 ++- cpp/purify/rm_kernel_integration.cc | 32 +++++++++++++++++++++++++++++ cpp/purify/rm_kernel_integration.h | 23 +++++++++++++++++++++ 3 files changed, 57 insertions(+), 1 deletion(-) create mode 100644 cpp/purify/rm_kernel_integration.cc create mode 100644 cpp/purify/rm_kernel_integration.h diff --git a/cpp/purify/CMakeLists.txt b/cpp/purify/CMakeLists.txt index 03a936a1c..f5255487b 100644 --- a/cpp/purify/CMakeLists.txt +++ b/cpp/purify/CMakeLists.txt @@ -29,12 +29,13 @@ set(HEADERS wkernel_integration.h wproj_operators.h rm_operators.h + rm_kernel_integration.h "${PROJECT_BINARY_DIR}/include/purify/config.h") set(SOURCES utilities.cc pfitsio.cc kernels.cc wproj_utilities.cc operators.cc uvfits.cc yaml-parser.cc read_measurements.cc distribute.cc integration.cc wide_field_utilities.cc wkernel_integration.cc - wproj_operators.cc rm_operators.cc) + wproj_operators.cc rm_operators.cc rm_kernel_integration.cc) if(TARGET casacore::ms) list(APPEND SOURCES casacore.cc) diff --git a/cpp/purify/rm_kernel_integration.cc b/cpp/purify/rm_kernel_integration.cc new file mode 100644 index 000000000..3f9d58ae4 --- /dev/null +++ b/cpp/purify/rm_kernel_integration.cc @@ -0,0 +1,32 @@ +#include "rm_kernel_integration.h" +#include + +namespace purify { + +namespace projection_kernels { +t_complex fourier_rm_kernel(const t_real x, const t_real dlambda2, const t_real u, + const t_real du) { + return std::exp(t_complex(0., -2 * constant::pi * u * x)) * + ((dlambda2 > 0) ? boost::math::sinc_pi(dlambda2 * x * constant::pi / du) : 1.); +} + +t_complex exact_rm_projection_integration(const t_real u, const t_real dlambda2, const t_real du, + const t_real oversample_ratio, + const std::function &ftkernelu, + const t_uint max_evaluations, + const t_real absolute_error, + const t_real relative_error, + const integration::method method, t_uint &evaluations) { + auto const func = [&](const Vector &x) -> t_complex { + evaluations++; + return ftkernelu(x(0)) * fourier_rm_kernel(x(0), dlambda2, u, du); + }; + Vector xmax = Vector::Zero(1); + xmax(0) = oversample_ratio / 2.; + const Vector xmin = -xmax; + return integration::integrate(xmin, xmax, func, integration::norm_type::paired, absolute_error, + relative_error, max_evaluations, method) / + (xmax(0) - xmin(0)); +} +} // namespace projection_kernels +} // namespace purify diff --git a/cpp/purify/rm_kernel_integration.h b/cpp/purify/rm_kernel_integration.h new file mode 100644 index 000000000..07bfc3598 --- /dev/null +++ b/cpp/purify/rm_kernel_integration.h @@ -0,0 +1,23 @@ +#ifndef PURIFY_RM_KERNEL_INTEGRATION_H +#define PURIFY_RM_KERNEL_INTEGRATION_H +#include "purify/config.h" +#include "purify/types.h" +#include "purify/logging.h" + +#include "purify/integration.h" +namespace purify { + +namespace projection_kernels { +//! fourier integrand with sinc function +t_complex fourier_rm_kernel(const t_real x, const t_real dlambda2, const t_real u, const t_real du); +//! integrate fourier integrand to get rm gridding kernel +t_complex exact_rm_projection_integration(const t_real u, const t_real dlambda2, const t_real du, + const t_real oversample_ratio, + const std::function &ftkernelu, + const t_uint max_evaluations, const t_real absolute_error, + const t_real relative_error, + const integration::method method, t_uint &evaluations); +} // namespace projection_kernels +// namespace projection_kernels +} // namespace purify +#endif From d0e9dbb015b1fca0528b02a1d93199793dc8b30a Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Tue, 23 Apr 2019 18:58:54 +0100 Subject: [PATCH 06/22] added gridding matrix that uses projection with rm kernels --- cpp/purify/rm_operators.cc | 60 ++++++++++++++++++++++++++++++++++++-- cpp/purify/rm_operators.h | 49 ++++++++++++++++--------------- 2 files changed, 84 insertions(+), 25 deletions(-) diff --git a/cpp/purify/rm_operators.cc b/cpp/purify/rm_operators.cc index 6a8d00cc4..35459e8fd 100644 --- a/cpp/purify/rm_operators.cc +++ b/cpp/purify/rm_operators.cc @@ -1,7 +1,11 @@ #include "purify/rm_operators.h" +#include "purify/rm_kernel_integration.h" namespace purify { -namespace details { +namespace rm_details { +t_real pixel_to_lambda2(const t_real cell, const t_uint imsize, const t_real oversample_ratio) { + return constant::pi / (oversample_ratio * cell * imsize); +} //! Construct gridding matrix Sparse init_gridding_matrix_1d(const Vector &u, const Vector &weights, const t_uint imsizex_, const t_real oversample_ratio, @@ -29,6 +33,58 @@ Sparse init_gridding_matrix_1d(const Vector &u, const Vector< } return interpolation_matrix; } +Sparse init_gridding_matrix_1d(const Vector &u, const Vector &widths, + const Vector &weights, const t_uint imsizex_, + const t_real cell, const t_real oversample_ratio, + const std::function ftkernelu, + const t_uint Ju, const t_uint J_max, + const t_real rel_error, const t_real abs_error) { + const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio); + const t_int rows = u.size(); + const t_int cols = ftsizeu_; + + const t_real du = rm_details::pixel_to_lambda2(cell, imsizex_, oversample_ratio); + + // count gridding coefficients with variable support size + Vector total_coeffs = Vector::Zero(widths.size()); + for (int i = 0; i < widths.size(); i++) { + const t_int Ju_max = std::min(std::floor(Ju + widths(i) / du), J_max); + total_coeffs(i) = Ju_max; + } + const t_real num_of_coeffs = total_coeffs.array().cast().sum(); + PURIFY_HIGH_LOG("Using {} rows (coefficients per a row {}), and memory of {} MegaBytes", + total_coeffs.size(), total_coeffs.array().cast().mean(), + 16. * num_of_coeffs / std::pow(10., 6)); + Sparse interpolation_matrix(rows, cols); + try { + interpolation_matrix.reserve(total_coeffs); + } catch (std::bad_alloc e) { + throw std::runtime_error( + "Not enough memory for coefficients, choose upper limit on support size Jw."); + } + interpolation_matrix.reserve(Vector::Constant(rows, Ju)); + + const t_complex I(0, 1); + const t_uint max_evaluations = 1e9; + // If I collapse this for loop there is a crash when using MPI... Sparse<>::insert() doesn't work + // right +#pragma omp parallel for + for (t_int m = 0; m < rows; ++m) { + const t_int ju_max = std::min(std::floor(Ju + widths(m) / du), J_max); + for (t_int ju = 1; ju < ju_max + 1; ++ju) { + const t_real k_u = std::floor(u(m) - ju_max * 0.5); + const t_uint index = utilities::mod(k_u + ju, ftsizeu_); + t_uint evaluations = 0; + interpolation_matrix.insert(m, index) = + std::exp(-2 * constant::pi * I * (k_u + ju) * 0.5) * + projection_kernels::exact_rm_projection_integration( + u(m) - (k_u + ju), widths(m), du, oversample_ratio, ftkernelu, max_evaluations, + abs_error, rel_error, integration::method::p, evaluations) * + weights(m); + } + } + return interpolation_matrix; +} Image init_correction1d(const t_real oversample_ratio, const t_uint imsizex_, const std::function ftkernelu) { @@ -42,5 +98,5 @@ Image init_correction1d(const t_real oversample_ratio, const t_uint i })) * t_complex(1., 0.); } -} // namespace details +} // namespace rm_details } // namespace purify diff --git a/cpp/purify/rm_operators.h b/cpp/purify/rm_operators.h index d051ce457..232b54bc0 100644 --- a/cpp/purify/rm_operators.h +++ b/cpp/purify/rm_operators.h @@ -8,19 +8,27 @@ namespace purify { -namespace details { +namespace rm_details { //! Construct gridding matrix Sparse init_gridding_matrix_1d(const Vector &u, const Vector &weights, const t_uint imsizex_, const t_real oversample_ratio, const std::function kernelu, const t_uint Ju); +//! Construct gridding matrix with fde +Sparse init_gridding_matrix_1d(const Vector &u, const Vector &widths, + const Vector &weights, const t_uint imsizex_, + const t_real cell, const t_real oversample_ratio, + const std::function ftkernelu, + const t_uint Ju, const t_uint J_max, + const t_real rel_error, const t_real abs_error); //! Given the Fourier transform of a gridding kernel, creates the scaling image for gridding //! correction. Image init_correction1d(const t_real oversample_ratio, const t_uint imsizex_, const std::function ftkernelu); - +//! gives size of lambda^2 pixel given faraday cell size (in rad/m^2) and image size +t_real pixel_to_lambda2(const t_real cell, const t_uint imsize, const t_real oversample_ratio); } // namespace details namespace operators { @@ -31,7 +39,7 @@ std::tuple, sopt::OperatorFunction> init_gridding_m ARGS &&... args) { const std::shared_ptr> interpolation_matrix = std::make_shared>( - details::init_gridding_matrix_1d(std::forward(args)...)); + rm_details::init_gridding_matrix_1d(std::forward(args)...)); const std::shared_ptr> adjoint = std::make_shared>(interpolation_matrix->adjoint()); @@ -135,7 +143,7 @@ std::tuple, sopt::OperatorFunction> base_padding_an sopt::OperatorFunction directFFT, indirectFFT; const Image S = - purify::details::init_correction1d(oversample_ratio, imsizex, ftkernelu); + purify::rm_details::init_correction1d(oversample_ratio, imsizex, ftkernelu); PURIFY_LOW_LOG("Building Measurement Operator: WGFZDB"); PURIFY_LOW_LOG( "Constructing Zero Padding " @@ -162,9 +170,11 @@ std::tuple, sopt::OperatorFunction> base_padding_an template std::tuple, sopt::OperatorFunction> base_degrid_operator_1d( - const Vector &u, const Vector &weights, const t_uint imsizex, - const t_real oversample_ratio = 2, const kernels::kernel kernel = kernels::kernel::kb, - const t_uint Ju = 4, const fftw_plan ft_plan = fftw_plan::measure) { + const Vector &u, const Vector &widths, const Vector &weights, + const t_uint imsizex, const t_real cell, const t_real oversample_ratio = 2, + const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4, + const t_uint J_max = 30, const t_real abs_error = 1e-6, const t_real rel_error = 1e-6, + const fftw_plan ft_plan = fftw_plan::measure) { std::function kernelu, ftkernelu; std::tie(ftkernelu, kernelu) = purify::create_radial_ftkernel(kernel, Ju, oversample_ratio); sopt::OperatorFunction directFZ, indirectFZ; @@ -174,7 +184,8 @@ std::tuple, sopt::OperatorFunction> base_degrid_ope PURIFY_LOW_LOG("Constructing Weighting and Gridding Operators: WG"); PURIFY_MEDIUM_LOG("Number of channels: {}", u.size()); std::tie(directG, indirectG) = purify::operators::init_gridding_matrix_1d( - u, weights, imsizex, oversample_ratio, kernelu, Ju); + u, widths, weights, imsizex, cell, oversample_ratio, kernelu, Ju, J_max, abs_error, + rel_error); auto direct = sopt::chained_operators(directG, directFZ); auto indirect = sopt::chained_operators(indirectFZ, indirectG); PURIFY_LOW_LOG("Finished consturction of RM Φ."); @@ -188,29 +199,21 @@ namespace measurementoperator { //! Returns linear transform that is the standard degridding operator template std::shared_ptr> init_degrid_operator_1d( - const Vector &u, const Vector &weights, const t_uint imsizex, - const t_real oversample_ratio = 2, const kernels::kernel kernel = kernels::kernel::kb, - const t_uint Ju = 4) { + const Vector &u, const Vector &widths, const Vector &weights, + const t_uint imsizex, const t_real cell, const t_real oversample_ratio = 2, + const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4, + const t_uint J_max = 30, const t_real abs_error = 1e-6, const t_real rel_error = 1e-6) { + const t_real du = rm_details::pixel_to_lambda2(cell, imsizex, oversample_ratio); const operators::fftw_plan ft_plan = operators::fftw_plan::measure; std::array N = {0, 1, static_cast(imsizex)}; std::array M = {0, 1, static_cast(u.size())}; sopt::OperatorFunction directDegrid, indirectDegrid; std::tie(directDegrid, indirectDegrid) = purify::operators::base_degrid_operator_1d( - u, weights, imsizex, oversample_ratio, kernel, Ju, ft_plan); + -u / du, widths, weights, imsizex, cell, oversample_ratio, kernel, Ju, J_max, abs_error, + rel_error, ft_plan); return std::make_shared>(directDegrid, M, indirectDegrid, N); } -template -std::shared_ptr> init_degrid_operator_1d( - const utilities::vis_params &uv_vis_input, const t_uint imsizex, - const t_real oversample_ratio = 2, const kernels::kernel kernel = kernels::kernel::kb, - const t_uint Ju = 4) { - auto uv_vis = uv_vis_input; - // need to find a way to set RM resolution, i.e. Faraday pixel size - return init_degrid_operator_1d(uv_vis.u, uv_vis.weights, imsizex, oversample_ratio, kernel, - Ju); -} - } // namespace measurementoperator }; // namespace purify #endif From 2899ad66ba96bb07f5b0468083c53b390b7773b9 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Tue, 23 Apr 2019 18:59:18 +0100 Subject: [PATCH 07/22] added example that plots rm kernel --- cpp/example/CMakeLists.txt | 1 + cpp/example/plot_rmkernel.cc | 96 ++++++++++++++++++++++++++++++++++++ 2 files changed, 97 insertions(+) create mode 100644 cpp/example/plot_rmkernel.cc diff --git a/cpp/example/CMakeLists.txt b/cpp/example/CMakeLists.txt index 51914e686..ad1fb0389 100644 --- a/cpp/example/CMakeLists.txt +++ b/cpp/example/CMakeLists.txt @@ -9,6 +9,7 @@ add_example(compare_wprojection LIBRARIES libpurify NOTEST) add_example(generate_vis_data LIBRARIES libpurify NOTEST) add_example(plot_wkernel LIBRARIES libpurify NOTEST) +add_example(plot_rmkernel LIBRARIES libpurify NOTEST) add_example(wavelet_decomposition LIBRARIES libpurify NOTEST) add_example(histogram_equalisation LIBRARIES libpurify NOTEST) diff --git a/cpp/example/plot_rmkernel.cc b/cpp/example/plot_rmkernel.cc new file mode 100644 index 000000000..d236569ed --- /dev/null +++ b/cpp/example/plot_rmkernel.cc @@ -0,0 +1,96 @@ + +#include "purify/config.h" + +#include +#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/rm_kernel_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(args[1]) : "kernel"; +#define ARGS_MACRO(NAME, ARGN, VALUE, TYPE) \ + auto const NAME = \ + static_cast((nargs > ARGN) ? std::stod(static_cast(args[ARGN])) : VALUE); + + ARGS_MACRO(dl2, 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) + +#undef ARGS_MACRO + purify::logging::initialize(); + purify::logging::set_level("debug"); + t_uint const J = 4; + t_int const Jl = 30; + const t_real oversample_ratio = 2; + PURIFY_LOW_LOG("image size: {}", imsize); + PURIFY_LOW_LOG("cell size: {}", cell); + const t_real du = constant::pi / (oversample_ratio * cell * imsize); + 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 = std::abs(projection_kernels::exact_rm_projection_integration( + 0, 0, du, oversample_ratio, normu, 1e9, 0, 1e-8, integration::method::p, e)); + auto const ftkernelu = [=](const t_real l) -> t_real { + return kernels::ft_kaiser_bessel(l, J) / norm; + }; + 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 Jl_max = 100; + Vector kernel = Vector::Zero(Ju_max * Jl_max); + Vector evals = Vector::Zero(Ju_max * Jl_max); + Vector support = Vector::Zero(Jl_max); + PURIFY_LOW_LOG("dlambda^2: {}", dl2); + PURIFY_LOW_LOG("Kernel size: {} x {}", Ju_max, Jl_max); + PURIFY_LOW_LOG("FoV (rad/m^2): {}", imsize * cell); + PURIFY_LOW_LOG("du: {}", du); +#pragma omp parallel for collapse(2) + for (int i = 0; i < Jl_max; i++) { + for (int j = 0; j < Ju_max; j++) { + t_uint evaluations = 0; + t_uint revaluations = 0; + t_uint e; + if (j / upsample < (dl2 * i / Jl_max / du + J) / 2) { + kernel(j * Jl_max + i) = projection_kernels::exact_rm_projection_integration( + j / upsample, dl2 * i / Jl_max, du, oversample_ratio, ftkernelu, max_evaluations, + absolute_error, 0, integration::method::p, evaluations); + evals(j * Jl_max + i) = evaluations; +#pragma omp critical(print) + { + PURIFY_LOW_LOG("u : {}, dlambda^2 : {}", j / upsample, dl2 * i / Jl_max); + PURIFY_LOW_LOG("Kernel value: {}", kernel(j * Jl_max + i)); + PURIFY_LOW_LOG("Evaluations: {}", evaluations); + } + if (kernel(j * Jl_max + i) != kernel(j * Jl_max + i)) + throw std::runtime_error("Kernel value is a nan."); + total += evaluations; + } + } + } + PURIFY_LOW_LOG("Total Evaluations: {}", total); + PURIFY_LOW_LOG("FoV (rad/m^2): {}", imsize * cell); + PURIFY_LOW_LOG("du: {}", du); + pfitsio::write2d(kernel.real(), Jl_max, Ju_max, output_name + "_real.fits"); + pfitsio::write2d(kernel.imag(), Jl_max, Ju_max, output_name + "_imag.fits"); + pfitsio::write2d(kernel.cwiseAbs(), Jl_max, Ju_max, output_name + ".fits"); + pfitsio::write2d(evals, Jl_max, Ju_max, output_name + "_evals.fits"); +} From bef011ec6c10595d597e7bf0e5394a427f087ea5 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Tue, 23 Apr 2019 18:59:34 +0100 Subject: [PATCH 08/22] updated test for rm operator --- cpp/tests/rm_operator.cc | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/cpp/tests/rm_operator.cc b/cpp/tests/rm_operator.cc index 17dbb69e0..1ecafb552 100644 --- a/cpp/tests/rm_operator.cc +++ b/cpp/tests/rm_operator.cc @@ -68,6 +68,7 @@ TEST_CASE("1d FFT") { TEST_CASE("regression_degrid") { const t_int imsizex = 256; + const t_real cell = constant::pi / imsizex; Vector u = Vector::Random(10) * imsizex / 2; Vector input = Vector::Zero(imsizex); input(static_cast(imsizex * 0.5)) = 1.; @@ -83,7 +84,7 @@ TEST_CASE("regression_degrid") { const kernels::kernel kernel = kernels::kernel::kb; const std::shared_ptr>> measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsizex, oversample_ratio, kernel, Ju); + u, Vector::Ones(M), imsizex, cell, oversample_ratio, kernel, Ju); const Vector y_test = *measure_op * input; CAPTURE(y_test.array() / y.array()); const t_real max_test = y_test.cwiseAbs().mean(); @@ -94,7 +95,7 @@ TEST_CASE("regression_degrid") { SECTION("pswf") { const kernels::kernel kernel = kernels::kernel::pswf; const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsizex, oversample_ratio, kernel, 6); + u, Vector::Ones(M), imsizex, cell, oversample_ratio, kernel, 6); const Vector y_test = *measure_op * input; CAPTURE(y_test.array() / y.array()); const t_real max_test = y_test.cwiseAbs().mean(); @@ -105,7 +106,7 @@ TEST_CASE("regression_degrid") { SECTION("gauss") { const kernels::kernel kernel = kernels::kernel::gauss; const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsizex, oversample_ratio, kernel, Ju); + u, Vector::Ones(M), imsizex, cell, oversample_ratio, kernel, Ju); const Vector y_test = *measure_op * input; CAPTURE(y_test.array() / y.array()); const t_real max_test = y_test.cwiseAbs().mean(); @@ -124,8 +125,9 @@ TEST_CASE("flux units") { const kernels::kernel kernel = kernels::kernel::kb; for (auto& J : {4, 5, 6, 7, 8}) { for (auto& imsize : {128, 256, 512}) { + const t_real cell = constant::pi / imsize; const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsize, oversample_ratio, kernel, J); + u, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, J); Vector input = Vector::Zero(imsize); input(static_cast(imsize * 0.5)) = 1.; const Vector y_test = @@ -147,12 +149,14 @@ TEST_CASE("flux units") { const kernels::kernel kernel = kernels::kernel::pswf; for (auto& J : {4, 5, 6, 7, 8}) { for (auto& imsize : {128, 256, 512}) { + const t_real cell = constant::pi / imsize; if (J != 6) CHECK_THROWS(measurementoperator::init_degrid_operator_1d>( - u * imsize / 2, Vector::Ones(M), imsize, oversample_ratio, kernel, J)); + u * imsize / 2, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, + J)); else { const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsize, oversample_ratio, kernel, J); + u, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, J); Vector input = Vector::Zero(imsize); input(static_cast(imsize * 0.5)) = 1.; const Vector y_test = @@ -175,8 +179,9 @@ TEST_CASE("flux units") { const kernels::kernel kernel = kernels::kernel::gauss; for (auto& J : {4, 5, 6, 7, 8}) { for (auto& imsize : {128, 256, 512}) { + const t_real cell = constant::pi / imsize; const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsize, oversample_ratio, kernel, J); + u, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, J); Vector input = Vector::Zero(imsize); input(static_cast(imsize * 0.5)) = 1.; const Vector y_test = @@ -189,8 +194,7 @@ TEST_CASE("flux units") { CAPTURE(imsize) CHECK(y_test.isApprox(y, 1e-2)); const Vector psf = - (measure_op->adjoint() * y) * - std::sqrt(input.size() * oversample_ratio ) / M; + (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M; CHECK(std::real(psf(static_cast(imsize * 0.5))) == Approx(1.).margin(0.01)); } } @@ -199,8 +203,9 @@ TEST_CASE("flux units") { const kernels::kernel kernel = kernels::kernel::box; for (auto& J : {4, 5, 6, 7, 8}) { for (auto& imsize : {128, 256, 512}) { + const t_real cell = constant::pi / imsize; const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsize, oversample_ratio, kernel, J); + u, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, J); Vector input = Vector::Zero(imsize); input(static_cast(imsize * 0.5)) = 1.; const Vector y_test = @@ -231,10 +236,12 @@ TEST_CASE("normed operator") { const kernels::kernel kernel = kernels::kernel::kb; for (auto& J : {4, 5, 6, 7, 8}) { for (auto& imsize : {128, 256, 512}) { + const t_real cell = constant::pi / imsize; const Vector init = Vector::Ones(imsize); auto power_method_result = sopt::algorithm::normalise_operator>( measurementoperator::init_degrid_operator_1d>( - u * imsize / 2, Vector::Ones(M), imsize, oversample_ratio, kernel, J), + u * imsize / 2, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, + J), power_iters, power_tol, init); const t_real norm = std::get<0>(power_method_result); const std::shared_ptr>> measure_op = From 43f55c44d409569d69c9f88c139b7865413af4b6 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Tue, 23 Apr 2019 23:16:41 +0100 Subject: [PATCH 09/22] fixed normalisation --- cpp/purify/rm_kernel_integration.cc | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/cpp/purify/rm_kernel_integration.cc b/cpp/purify/rm_kernel_integration.cc index 3f9d58ae4..f1f3ee16e 100644 --- a/cpp/purify/rm_kernel_integration.cc +++ b/cpp/purify/rm_kernel_integration.cc @@ -13,8 +13,7 @@ t_complex fourier_rm_kernel(const t_real x, const t_real dlambda2, const t_real t_complex exact_rm_projection_integration(const t_real u, const t_real dlambda2, const t_real du, const t_real oversample_ratio, const std::function &ftkernelu, - const t_uint max_evaluations, - const t_real absolute_error, + const t_uint max_evaluations, const t_real absolute_error, const t_real relative_error, const integration::method method, t_uint &evaluations) { auto const func = [&](const Vector &x) -> t_complex { @@ -25,8 +24,7 @@ t_complex exact_rm_projection_integration(const t_real u, const t_real dlambda2, xmax(0) = oversample_ratio / 2.; const Vector xmin = -xmax; return integration::integrate(xmin, xmax, func, integration::norm_type::paired, absolute_error, - relative_error, max_evaluations, method) / - (xmax(0) - xmin(0)); + relative_error, max_evaluations, method); } } // namespace projection_kernels } // namespace purify From 6ed212fb1a0e37b1a95178f85bc5433c30ee190a Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Tue, 23 Apr 2019 23:17:03 +0100 Subject: [PATCH 10/22] formatting --- cpp/purify/rm_operators.cc | 2 -- cpp/purify/rm_operators.h | 11 +++++------ 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/cpp/purify/rm_operators.cc b/cpp/purify/rm_operators.cc index 35459e8fd..b15890ce0 100644 --- a/cpp/purify/rm_operators.cc +++ b/cpp/purify/rm_operators.cc @@ -66,8 +66,6 @@ Sparse init_gridding_matrix_1d(const Vector &u, const Vector< const t_complex I(0, 1); const t_uint max_evaluations = 1e9; - // If I collapse this for loop there is a crash when using MPI... Sparse<>::insert() doesn't work - // right #pragma omp parallel for for (t_int m = 0; m < rows; ++m) { const t_int ju_max = std::min(std::floor(Ju + widths(m) / du), J_max); diff --git a/cpp/purify/rm_operators.h b/cpp/purify/rm_operators.h index 232b54bc0..ab3fbb75d 100644 --- a/cpp/purify/rm_operators.h +++ b/cpp/purify/rm_operators.h @@ -29,7 +29,7 @@ Image init_correction1d(const t_real oversample_ratio, const t_uint i const std::function ftkernelu); //! gives size of lambda^2 pixel given faraday cell size (in rad/m^2) and image size t_real pixel_to_lambda2(const t_real cell, const t_uint imsize, const t_real oversample_ratio); -} // namespace details +} // namespace rm_details namespace operators { @@ -171,10 +171,9 @@ std::tuple, sopt::OperatorFunction> base_padding_an template std::tuple, sopt::OperatorFunction> base_degrid_operator_1d( const Vector &u, const Vector &widths, const Vector &weights, - const t_uint imsizex, const t_real cell, const t_real oversample_ratio = 2, - const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4, - const t_uint J_max = 30, const t_real abs_error = 1e-6, const t_real rel_error = 1e-6, - const fftw_plan ft_plan = fftw_plan::measure) { + const t_uint imsizex, const t_real cell, const t_real oversample_ratio, + const kernels::kernel kernel, const t_uint Ju, const t_uint J_max, const t_real abs_error, + const t_real rel_error, const fftw_plan ft_plan) { std::function kernelu, ftkernelu; std::tie(ftkernelu, kernelu) = purify::create_radial_ftkernel(kernel, Ju, oversample_ratio); sopt::OperatorFunction directFZ, indirectFZ; @@ -184,7 +183,7 @@ std::tuple, sopt::OperatorFunction> base_degrid_ope PURIFY_LOW_LOG("Constructing Weighting and Gridding Operators: WG"); PURIFY_MEDIUM_LOG("Number of channels: {}", u.size()); std::tie(directG, indirectG) = purify::operators::init_gridding_matrix_1d( - u, widths, weights, imsizex, cell, oversample_ratio, kernelu, Ju, J_max, abs_error, + u, widths, weights, imsizex, cell, oversample_ratio, ftkernelu, Ju, J_max, abs_error, rel_error); auto direct = sopt::chained_operators(directG, directFZ); auto indirect = sopt::chained_operators(indirectFZ, indirectG); From 8db7a176c60d0405c72260943d3d52677ce79764 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Tue, 23 Apr 2019 23:17:15 +0100 Subject: [PATCH 11/22] updated test --- cpp/tests/rm_operator.cc | 84 ++++++++-------------------------------- 1 file changed, 16 insertions(+), 68 deletions(-) diff --git a/cpp/tests/rm_operator.cc b/cpp/tests/rm_operator.cc index 1ecafb552..99520b366 100644 --- a/cpp/tests/rm_operator.cc +++ b/cpp/tests/rm_operator.cc @@ -69,7 +69,8 @@ TEST_CASE("1d FFT") { TEST_CASE("regression_degrid") { const t_int imsizex = 256; const t_real cell = constant::pi / imsizex; - Vector u = Vector::Random(10) * imsizex / 2; + const Vector u = Vector::Random(10) * imsizex / 2; + const Vector widths = Vector::Zero(u.size()); Vector input = Vector::Zero(imsizex); input(static_cast(imsizex * 0.5)) = 1.; const t_uint M = u.size(); @@ -84,7 +85,8 @@ TEST_CASE("regression_degrid") { const kernels::kernel kernel = kernels::kernel::kb; const std::shared_ptr>> measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsizex, cell, oversample_ratio, kernel, Ju); + u, widths, Vector::Ones(M), imsizex, cell, oversample_ratio, kernel, Ju, 30, + 1e-6, 1e-6); const Vector y_test = *measure_op * input; CAPTURE(y_test.array() / y.array()); const t_real max_test = y_test.cwiseAbs().mean(); @@ -95,7 +97,8 @@ TEST_CASE("regression_degrid") { SECTION("pswf") { const kernels::kernel kernel = kernels::kernel::pswf; const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsizex, cell, oversample_ratio, kernel, 6); + u, widths, Vector::Ones(M), imsizex, cell, oversample_ratio, kernel, 6, 30, 1e-6, + 1e-6); const Vector y_test = *measure_op * input; CAPTURE(y_test.array() / y.array()); const t_real max_test = y_test.cwiseAbs().mean(); @@ -103,22 +106,12 @@ TEST_CASE("regression_degrid") { CAPTURE(y); CHECK((y_test / max_test).isApprox((y), 1e-3)); } - SECTION("gauss") { - const kernels::kernel kernel = kernels::kernel::gauss; - const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsizex, cell, oversample_ratio, kernel, Ju); - const Vector y_test = *measure_op * input; - CAPTURE(y_test.array() / y.array()); - const t_real max_test = y_test.cwiseAbs().mean(); - CAPTURE(y_test / max_test); - CAPTURE(y); - CHECK((y_test / max_test).isApprox((y), 1e-4)); - } } TEST_CASE("flux units") { const t_real oversample_ratio = 2; Vector u = Vector::Random(10); + const Vector widths = Vector::Zero(u.size()); const t_uint M = u.size(); const Vector y = Vector::Ones(u.size()); SECTION("kb") { @@ -127,7 +120,8 @@ TEST_CASE("flux units") { for (auto& imsize : {128, 256, 512}) { const t_real cell = constant::pi / imsize; const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, J); + u * imsize / 2, widths, Vector::Ones(M), imsize, cell, oversample_ratio, + kernel, J, 30, 1e-6, 1e-6); Vector input = Vector::Zero(imsize); input(static_cast(imsize * 0.5)) = 1.; const Vector y_test = @@ -152,11 +146,12 @@ TEST_CASE("flux units") { const t_real cell = constant::pi / imsize; if (J != 6) CHECK_THROWS(measurementoperator::init_degrid_operator_1d>( - u * imsize / 2, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, - J)); + u * imsize / 2, widths, Vector::Ones(M), imsize, cell, oversample_ratio, + kernel, J, 30, 1e-6, 1e-6)); else { const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, J); + u * imsize / 2, widths, Vector::Ones(M), imsize, cell, oversample_ratio, + kernel, J); Vector input = Vector::Zero(imsize); input(static_cast(imsize * 0.5)) = 1.; const Vector y_test = @@ -175,54 +170,6 @@ TEST_CASE("flux units") { } } } - SECTION("gauss") { - const kernels::kernel kernel = kernels::kernel::gauss; - for (auto& J : {4, 5, 6, 7, 8}) { - for (auto& imsize : {128, 256, 512}) { - const t_real cell = constant::pi / imsize; - const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, J); - Vector input = Vector::Zero(imsize); - input(static_cast(imsize * 0.5)) = 1.; - const Vector y_test = - (*measure_op * input).eval() * - std::sqrt(static_cast(input.size()) * oversample_ratio); - CAPTURE(y_test.cwiseAbs().mean()); - CAPTURE(y); - CAPTURE(y_test); - CAPTURE(J); - CAPTURE(imsize) - CHECK(y_test.isApprox(y, 1e-2)); - const Vector psf = - (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M; - CHECK(std::real(psf(static_cast(imsize * 0.5))) == Approx(1.).margin(0.01)); - } - } - } - SECTION("box") { - const kernels::kernel kernel = kernels::kernel::box; - for (auto& J : {4, 5, 6, 7, 8}) { - for (auto& imsize : {128, 256, 512}) { - const t_real cell = constant::pi / imsize; - const auto measure_op = measurementoperator::init_degrid_operator_1d>( - u, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, J); - Vector input = Vector::Zero(imsize); - input(static_cast(imsize * 0.5)) = 1.; - const Vector y_test = - (*measure_op * input).eval() * - std::sqrt(static_cast(input.size()) * oversample_ratio); - CAPTURE(y_test.cwiseAbs().mean()); - CAPTURE(y); - CAPTURE(y_test); - CAPTURE(J); - CAPTURE(imsize) - CHECK(y_test.isApprox(y, 1e-3)); - const Vector psf = - (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M; - CHECK(std::real(psf(static_cast(imsize * 0.5))) == Approx(1.).margin(0.001)); - } - } - } } TEST_CASE("normed operator") { @@ -230,6 +177,7 @@ TEST_CASE("normed operator") { const t_int power_iters = 1000; const t_real power_tol = 1e-4; Vector u = Vector::Random(10); + const Vector widths = Vector::Zero(u.size()); const t_uint M = u.size(); const Vector y = Vector::Ones(u.size()); SECTION("kb") { @@ -240,8 +188,8 @@ TEST_CASE("normed operator") { const Vector init = Vector::Ones(imsize); auto power_method_result = sopt::algorithm::normalise_operator>( measurementoperator::init_degrid_operator_1d>( - u * imsize / 2, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, - J), + u * imsize / 2, widths, Vector::Ones(M), imsize, cell, oversample_ratio, + kernel, J, 30, 1e-6, 1e-6), power_iters, power_tol, init); const t_real norm = std::get<0>(power_method_result); const std::shared_ptr>> measure_op = From 89bc56166cb0dc901ae20fe13bce30293e6f16c9 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Fri, 3 May 2019 16:34:32 +0100 Subject: [PATCH 12/22] added print out information for rm operator --- cpp/purify/rm_operators.cc | 1 + cpp/purify/rm_operators.h | 7 +++++++ 2 files changed, 8 insertions(+) diff --git a/cpp/purify/rm_operators.cc b/cpp/purify/rm_operators.cc index b15890ce0..3260bc6d3 100644 --- a/cpp/purify/rm_operators.cc +++ b/cpp/purify/rm_operators.cc @@ -52,6 +52,7 @@ Sparse init_gridding_matrix_1d(const Vector &u, const Vector< total_coeffs(i) = Ju_max; } const t_real num_of_coeffs = total_coeffs.array().cast().sum(); + PURIFY_HIGH_LOG("Max channel width in wavelengths squared: {}", widths.maxCoeff()); PURIFY_HIGH_LOG("Using {} rows (coefficients per a row {}), and memory of {} MegaBytes", total_coeffs.size(), total_coeffs.array().cast().mean(), 16. * num_of_coeffs / std::pow(10., 6)); diff --git a/cpp/purify/rm_operators.h b/cpp/purify/rm_operators.h index ab3fbb75d..6df6bb829 100644 --- a/cpp/purify/rm_operators.h +++ b/cpp/purify/rm_operators.h @@ -202,6 +202,13 @@ std::shared_ptr> init_degrid_operator_1d( const t_uint imsizex, const t_real cell, const t_real oversample_ratio = 2, const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4, const t_uint J_max = 30, const t_real abs_error = 1e-6, const t_real rel_error = 1e-6) { + PURIFY_LOW_LOG("Cell size in Faraday Depth is {} rad/m^2", cell); + const t_real min_sensitivity = 1.895 / widths.maxCoeff(); + const t_real max_sensitivity = 1.895 / widths.minCoeff(); + const t_real imaging_sensitivity = imsizex * cell / 2; + PURIFY_LOW_LOG("The max range of sensitivity is +/- {} rad/m^2", max_sensitivity); + PURIFY_LOW_LOG("The min range of sensitivity is +/- {} rad/m^2", min_sensitivity); + PURIFY_LOW_LOG("The imaging range is +/- {} rad/m^2", imsizex * cell / 2); const t_real du = rm_details::pixel_to_lambda2(cell, imsizex, oversample_ratio); const operators::fftw_plan ft_plan = operators::fftw_plan::measure; std::array N = {0, 1, static_cast(imsizex)}; From e1d052d04d00c517cdfd23e836a71ece7c700aca Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Fri, 3 May 2019 16:35:28 +0100 Subject: [PATCH 13/22] added rm reconstruction example with and without channel correction --- cpp/example/rm_gridding_example.cc | 131 +++++++++++++++++++++++++++++ 1 file changed, 131 insertions(+) create mode 100644 cpp/example/rm_gridding_example.cc diff --git a/cpp/example/rm_gridding_example.cc b/cpp/example/rm_gridding_example.cc new file mode 100644 index 000000000..1942dfb31 --- /dev/null +++ b/cpp/example/rm_gridding_example.cc @@ -0,0 +1,131 @@ +#include "purify/config.h" +#include "purify/types.h" +#include "purify/logging.h" +#include "purify/pfitsio.h" +#include "purify/rm_operators.h" +#include "purify/utilities.h" +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace purify; + +namespace { +void run_admm(const Vector &y, + const std::shared_ptr>> &measure_op, + const sopt::LinearTransform> &Psi, const t_real epsilon, + const std::string &prefix, const t_real oversample_ratio, const t_real channels, + const t_real imsizex, const t_real measure_op_norm) { + const Vector dmap = (measure_op->adjoint() * y) / channels * measure_op_norm * + measure_op_norm * oversample_ratio * imsizex; // in jy/beam + pfitsio::write2d(dmap.real(), prefix + "_rm_dmap_real.fits"); + pfitsio::write2d(dmap.imag(), prefix + "_rm_dmap_imag.fits"); + + auto const padmm = + sopt::algorithm::ImagingProximalADMM(y) + .itermax(10000) + .gamma((Psi.adjoint() * (measure_op->adjoint() * y).eval()).cwiseAbs().maxCoeff() * 1e-3) + .relative_variation(1e-6) + .l2ball_proximal_epsilon(epsilon) + .tight_frame(false) + .l1_proximal_tolerance(1e-3) + .l1_proximal_nu(1.) + .l1_proximal_itermax(50) + .l1_proximal_positivity_constraint(false) + .l1_proximal_real_constraint(false) + .residual_convergence(epsilon) + .lagrange_update_scale(0.9) + .nu(1e0) + .Psi(Psi) + .Phi(*measure_op); + + auto const diagnostic = padmm(); + Image image = Image::Map(diagnostic.x.data(), 1, imsizex); + pfitsio::write2d(image.real(), prefix + "_sol_real.fits"); + pfitsio::write2d(image.imag(), prefix + "_sol_imag.fits"); + Vector residuals = measure_op->adjoint() * (y - ((*measure_op) * diagnostic.x)) / + channels * measure_op_norm * measure_op_norm * oversample_ratio * + imsizex; + Image residual_image = Image::Map(residuals.data(), 1, imsizex); + pfitsio::write2d(residual_image.real(), prefix + "_res_real.fits"); + pfitsio::write2d(residual_image.imag(), prefix + "_res_imag.fits"); +} +} // namespace + +int main(int nargs, char const **args) { + purify::logging::initialize(); + purify::logging::set_level("debug"); + sopt::logging::initialize(); + sopt::logging::set_level("debug"); + + const t_int channels = 1100; + const t_real start_c = 700e6; + const t_real end_c = 1800e6; + const t_real snr = 30; + auto const oversample_ratio = 2; + auto const power_iters = 100; + auto const power_tol = 1e-4; + auto const Ju = 4; + auto const imsizex = 8192; + + auto const kernel = "kb"; + const t_real dnu = (end_c - start_c) / channels; // Hz + const t_int Jl_max = 100; + const Vector freq = Vector::LinSpaced(channels, start_c, end_c); + const Vector lambda2 = ((constant::c / (freq.array() - dnu)).square() + + (constant::c / (freq.array() + dnu)).square()) / + 2; + const Vector widths = ((constant::c / (freq.array() - dnu)).square() - + (constant::c / (freq.array() + dnu)).square()); + + const t_real cell = constant::pi / (2 * lambda2.maxCoeff()); // rad/m^2 + const auto measure_op_stuff = sopt::algorithm::normalise_operator>( + measurementoperator::init_degrid_operator_1d>( + lambda2, widths, Vector::Ones(lambda2.size()), imsizex, cell, oversample_ratio, + kernels::kernel_from_string.at(kernel), Ju, Jl_max, 1e-6, 1e-6), + power_iters, power_tol, Vector::Random(imsizex).eval()); + const auto measure_op = std::get<2>(measure_op_stuff); + const t_real measure_op_norm = std::get<0>(measure_op_stuff); + const t_real scale = std::sqrt(imsizex * oversample_ratio); + + const auto measure_op_no_correction_stuff = + sopt::algorithm::normalise_operator>( + measurementoperator::init_degrid_operator_1d>( + lambda2, Vector::Zero(lambda2.size()), + Vector::Ones(lambda2.size()), imsizex, cell, oversample_ratio, + kernels::kernel_from_string.at(kernel), Ju, Jl_max, 1e-6, 1e-6), + power_iters, power_tol, Vector::Random(imsizex).eval()); + const auto measure_op_no_correction = std::get<2>(measure_op_no_correction_stuff); + const t_real measure_op_norm_no_correction = std::get<0>(measure_op_no_correction_stuff); + + Vector input = Vector::Zero(imsizex); + for (t_int i = 0; i < 15; i++) + input((i + 1) * 512 - 1) = + std::exp(-2. * t_complex(0., 1.) * constant::pi * static_cast(i) / 8.); + + Vector y = (*measure_op * input); + t_real const sigma = utilities::SNR_to_standard_deviation(y, snr); + + // adding noise to visibilities + y = utilities::add_noise(y, 0., sigma); + auto const epsilon = utilities::calculate_l2_radius(y.size(), sigma); + + sopt::wavelets::SARA const sara{std::make_tuple("Dirac", 3u)}; + //, std::make_tuple("DB1", 3u), std::make_tuple("DB2", 3u), + // std::make_tuple("DB3", 3u), std::make_tuple("DB4", 3u), std::make_tuple("DB5", 3u), + // std::make_tuple("DB6", 3u), std::make_tuple("DB7", 3u), + // std::make_tuple("DB8", 3u)}; + + auto const Psi = sopt::linear_transform(sara, imsizex, 1); + + run_admm(y, measure_op, Psi, epsilon, "corrected", oversample_ratio, channels, imsizex, + measure_op_norm); + run_admm(y, measure_op_no_correction, Psi, epsilon, "not_corrected", oversample_ratio, channels, + imsizex, measure_op_norm_no_correction); +} From 2ba20348e7fe18c57679fc9eaf5df3a69bf66381 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Fri, 3 May 2019 16:38:07 +0100 Subject: [PATCH 14/22] added an example that will compare rm kernel correction against sinc function --- cpp/example/image_rm_window.cc | 65 ++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 cpp/example/image_rm_window.cc diff --git a/cpp/example/image_rm_window.cc b/cpp/example/image_rm_window.cc new file mode 100644 index 000000000..d0ab3c49e --- /dev/null +++ b/cpp/example/image_rm_window.cc @@ -0,0 +1,65 @@ +#include "purify/config.h" +#include "purify/types.h" +#include "purify/directories.h" +#include "purify/logging.h" +#include "purify/operators.h" +#include "purify/pfitsio.h" +#include "purify/rm_operators.h" +#include "purify/utilities.h" +#include + +#include + +using namespace purify; +using namespace purify::notinstalled; + +int main(int nargs, char const **args) { +#define ARGS_MACRO(NAME, ARGN, VALUE, TYPE) \ + auto const NAME = \ + static_cast((nargs > ARGN) ? std::stod(static_cast(args[ARGN])) : VALUE); + + ARGS_MACRO(width, 1, 10, t_real) + ARGS_MACRO(imsize, 2, 256, t_uint) + ARGS_MACRO(cell, 3, 2400, t_real) + ARGS_MACRO(Jl_max, 4, 0, t_uint) +#undef ARGS_MACRO + purify::logging::initialize(); + purify::logging::set_level("debug"); + // Gridding example + auto const oversample_ratio = 2; + auto const power_iters = 100; + auto const power_tol = 1e-5; + auto const Ju = 4; + t_int const imsizex = imsize; + const std::string suffix = "_" + std::to_string(static_cast(width)) + "_" + + std::to_string(static_cast(imsize)) + "_" + + std::to_string(static_cast(cell)); + + auto const kernel = "kb"; + + t_uint const number_of_pixels = imsizex; + t_uint const number_of_vis = std::floor(number_of_pixels * 2); + // Generating random uv(w) coverage + t_real const sigma_m = constant::pi / 3; + + auto header = + pfitsio::header_params("", " ", 1, 0, 0, stokes::I, cell, cell, 0, 1, 0, 0, 0, 0, 0); + const auto measure_opw = measurementoperator::init_degrid_operator_1d>( + Vector::Zero(1), Vector::Constant(1, width), Vector::Ones(1), + imsizex, cell, oversample_ratio, kernels::kernel_from_string.at(kernel), Ju, Jl_max, 1e-6, + 1e-6); + Vector const measurements = + (measure_opw->adjoint() * Vector::Constant(1, 1)).conjugate(); + t_uint max_pos = 0; + measurements.array().real().maxCoeff(&max_pos); + Image out = measurements * std::sqrt(imsizex * oversample_ratio); + header.fits_name = "rm_window_real" + suffix + ".fits"; + pfitsio::write2d(out.real(), header); + header.fits_name = "rm_window_imag" + suffix + ".fits"; + pfitsio::write2d(out.imag(), header); + Vector true_window = Vector::Zero(imsize); + for (t_int i = 0; i < imsize; i++) + true_window(i) = boost::math::sinc_pi(width * std::abs(i - (imsize) / 2.) * cell); + header.fits_name = "true_window_real" + suffix + ".fits"; + pfitsio::write2d(true_window.real(), header); +} From 53099ed158d6a7c8d66f53293da44f6ab2eb1867 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Fri, 3 May 2019 16:39:03 +0100 Subject: [PATCH 15/22] updated cmake to compile examples --- cpp/example/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cpp/example/CMakeLists.txt b/cpp/example/CMakeLists.txt index ad1fb0389..a32b53a5c 100644 --- a/cpp/example/CMakeLists.txt +++ b/cpp/example/CMakeLists.txt @@ -4,7 +4,9 @@ include_directories(SYSTEM ${Sopt_INCLUDE_DIRS}) include_directories("${PROJECT_SOURCE_DIR}/cpp") add_example(gridding LIBRARIES libpurify NOTEST) +add_example(rm_gridding_example LIBRARIES libpurify NOTEST) add_example(image_wproj_chirp LIBRARIES libpurify NOTEST) +add_example(image_rm_window LIBRARIES libpurify NOTEST) add_example(compare_wprojection LIBRARIES libpurify NOTEST) add_example(generate_vis_data LIBRARIES libpurify NOTEST) From ff67efaa49bbbf6804341be0ee67aee367ac07d9 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Fri, 3 May 2019 16:40:12 +0100 Subject: [PATCH 16/22] updated support size in example that generates rm kernel plots --- cpp/example/plot_rmkernel.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/example/plot_rmkernel.cc b/cpp/example/plot_rmkernel.cc index d236569ed..7ae9968b8 100644 --- a/cpp/example/plot_rmkernel.cc +++ b/cpp/example/plot_rmkernel.cc @@ -52,7 +52,7 @@ int main(int nargs, char const **args) { }; t_uint total = 0; t_uint rtotal = 0; - t_int const Ju = 20; + t_int const Ju = 40; t_real const upsample = 5; t_int const Ju_max = std::floor(Ju * upsample); const t_int Jl_max = 100; From 891a26005e1b09302ddc479be74ea491017f6c58 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Wed, 8 May 2019 00:18:00 +0100 Subject: [PATCH 17/22] fixed lambda calculation --- cpp/example/rm_gridding_example.cc | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/cpp/example/rm_gridding_example.cc b/cpp/example/rm_gridding_example.cc index 1942dfb31..9f1d4d938 100644 --- a/cpp/example/rm_gridding_example.cc +++ b/cpp/example/rm_gridding_example.cc @@ -29,9 +29,9 @@ void run_admm(const Vector &y, auto const padmm = sopt::algorithm::ImagingProximalADMM(y) - .itermax(10000) + .itermax(100000) .gamma((Psi.adjoint() * (measure_op->adjoint() * y).eval()).cwiseAbs().maxCoeff() * 1e-3) - .relative_variation(1e-6) + .relative_variation(1e-7) .l2ball_proximal_epsilon(epsilon) .tight_frame(false) .l1_proximal_tolerance(1e-3) @@ -72,17 +72,17 @@ int main(int nargs, char const **args) { auto const power_iters = 100; auto const power_tol = 1e-4; auto const Ju = 4; - auto const imsizex = 8192; + auto const imsizex = 8192 * 2; auto const kernel = "kb"; const t_real dnu = (end_c - start_c) / channels; // Hz const t_int Jl_max = 100; const Vector freq = Vector::LinSpaced(channels, start_c, end_c); - const Vector lambda2 = ((constant::c / (freq.array() - dnu)).square() + - (constant::c / (freq.array() + dnu)).square()) / + const Vector lambda2 = ((constant::c / (freq.array() - dnu / 2.)).square() + + (constant::c / (freq.array() + dnu / 2.)).square()) / 2; - const Vector widths = ((constant::c / (freq.array() - dnu)).square() - - (constant::c / (freq.array() + dnu)).square()); + const Vector widths = ((constant::c / (freq.array() - dnu / 2.)).square() - + (constant::c / (freq.array() + dnu / 2.)).square()); const t_real cell = constant::pi / (2 * lambda2.maxCoeff()); // rad/m^2 const auto measure_op_stuff = sopt::algorithm::normalise_operator>( @@ -106,8 +106,10 @@ int main(int nargs, char const **args) { Vector input = Vector::Zero(imsizex); for (t_int i = 0; i < 15; i++) - input((i + 1) * 512 - 1) = + input((i + 1) * 1024 - 1) = std::exp(-2. * t_complex(0., 1.) * constant::pi * static_cast(i) / 8.); + pfitsio::write2d(input.real(), "rm_input_real.fits"); + pfitsio::write2d(input.imag(), "rm_input_imag.fits"); Vector y = (*measure_op * input); t_real const sigma = utilities::SNR_to_standard_deviation(y, snr); From 2e8f28d6adcc4762a033eac77f4f10789f7ba7d0 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Wed, 8 May 2019 00:23:31 +0100 Subject: [PATCH 18/22] changed support of max kernel plot back --- cpp/example/plot_rmkernel.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/example/plot_rmkernel.cc b/cpp/example/plot_rmkernel.cc index 7ae9968b8..d236569ed 100644 --- a/cpp/example/plot_rmkernel.cc +++ b/cpp/example/plot_rmkernel.cc @@ -52,7 +52,7 @@ int main(int nargs, char const **args) { }; t_uint total = 0; t_uint rtotal = 0; - t_int const Ju = 40; + t_int const Ju = 20; t_real const upsample = 5; t_int const Ju_max = std::floor(Ju * upsample); const t_int Jl_max = 100; From fa45cb7dd3ecc50ebc5e3554047f079baf74a241 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Wed, 5 Jun 2019 18:34:59 +0100 Subject: [PATCH 19/22] added more information about RM sensitivity range --- cpp/purify/rm_operators.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/cpp/purify/rm_operators.h b/cpp/purify/rm_operators.h index 6df6bb829..dbea300c6 100644 --- a/cpp/purify/rm_operators.h +++ b/cpp/purify/rm_operators.h @@ -205,9 +205,15 @@ std::shared_ptr> init_degrid_operator_1d( PURIFY_LOW_LOG("Cell size in Faraday Depth is {} rad/m^2", cell); const t_real min_sensitivity = 1.895 / widths.maxCoeff(); const t_real max_sensitivity = 1.895 / widths.minCoeff(); + + const t_real min_null = constant::pi / widths.maxCoeff(); + const t_real max_null = constant::pi / widths.minCoeff(); + const t_real imaging_sensitivity = imsizex * cell / 2; PURIFY_LOW_LOG("The max range of sensitivity is +/- {} rad/m^2", max_sensitivity); PURIFY_LOW_LOG("The min range of sensitivity is +/- {} rad/m^2", min_sensitivity); + PURIFY_LOW_LOG("The last sensitivity null is +/- {} rad/m^2", max_null); + PURIFY_LOW_LOG("The first sensitivity null is +/- {} rad/m^2", min_null); PURIFY_LOW_LOG("The imaging range is +/- {} rad/m^2", imsizex * cell / 2); const t_real du = rm_details::pixel_to_lambda2(cell, imsizex, oversample_ratio); const operators::fftw_plan ft_plan = operators::fftw_plan::measure; From 8911119ff11fa0833d95a337b8145e093e2e18f2 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Wed, 5 Jun 2019 18:35:43 +0100 Subject: [PATCH 20/22] added MWA and POSSUM coverages to RM example --- cpp/example/rm_gridding_example.cc | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/cpp/example/rm_gridding_example.cc b/cpp/example/rm_gridding_example.cc index 9f1d4d938..b41072864 100644 --- a/cpp/example/rm_gridding_example.cc +++ b/cpp/example/rm_gridding_example.cc @@ -63,11 +63,24 @@ int main(int nargs, char const **args) { purify::logging::set_level("debug"); sopt::logging::initialize(); sopt::logging::set_level("debug"); - +/* POSSUM const t_int channels = 1100; const t_real start_c = 700e6; const t_real end_c = 1800e6; const t_real snr = 30; + const t_real dnu = 1e6; // Hz + auto const oversample_ratio = 2; + auto const power_iters = 100; + auto const power_tol = 1e-4; + auto const Ju = 4; + auto const imsizex = 8192 * 2; +*/ + + const t_int channels = 768; + const t_real start_c = 200.32e6; + const t_real end_c = 231.04e6; + const t_real snr = 30; + const t_real dnu = 40e3; // Hz auto const oversample_ratio = 2; auto const power_iters = 100; auto const power_tol = 1e-4; @@ -75,7 +88,6 @@ int main(int nargs, char const **args) { auto const imsizex = 8192 * 2; auto const kernel = "kb"; - const t_real dnu = (end_c - start_c) / channels; // Hz const t_int Jl_max = 100; const Vector freq = Vector::LinSpaced(channels, start_c, end_c); const Vector lambda2 = ((constant::c / (freq.array() - dnu / 2.)).square() + From 5f5189a827bda43b72f0cd8a65ba8ea37abf4d94 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Fri, 19 Jul 2019 15:03:31 +0100 Subject: [PATCH 21/22] linting --- cpp/example/rm_gridding_example.cc | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/cpp/example/rm_gridding_example.cc b/cpp/example/rm_gridding_example.cc index b41072864..0b5c89383 100644 --- a/cpp/example/rm_gridding_example.cc +++ b/cpp/example/rm_gridding_example.cc @@ -63,18 +63,18 @@ int main(int nargs, char const **args) { purify::logging::set_level("debug"); sopt::logging::initialize(); sopt::logging::set_level("debug"); -/* POSSUM - const t_int channels = 1100; - const t_real start_c = 700e6; - const t_real end_c = 1800e6; - const t_real snr = 30; - const t_real dnu = 1e6; // Hz - auto const oversample_ratio = 2; - auto const power_iters = 100; - auto const power_tol = 1e-4; - auto const Ju = 4; - auto const imsizex = 8192 * 2; -*/ + /* POSSUM + const t_int channels = 1100; + const t_real start_c = 700e6; + const t_real end_c = 1800e6; + const t_real snr = 30; + const t_real dnu = 1e6; // Hz + auto const oversample_ratio = 2; + auto const power_iters = 100; + auto const power_tol = 1e-4; + auto const Ju = 4; + auto const imsizex = 8192 * 2; + */ const t_int channels = 768; const t_real start_c = 200.32e6; From 34445ee1707c787cc1ff7152fede9f0ec3ba8d15 Mon Sep 17 00:00:00 2001 From: Luke Pratley Date: Fri, 4 Oct 2019 00:30:32 +0200 Subject: [PATCH 22/22] added rm operator with no channel averaging built in, with tests --- cpp/purify/rm_operators.h | 38 +++++++++++++++++++++++++++++++++++ cpp/tests/rm_operator.cc | 42 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 80 insertions(+) diff --git a/cpp/purify/rm_operators.h b/cpp/purify/rm_operators.h index dbea300c6..1115b2338 100644 --- a/cpp/purify/rm_operators.h +++ b/cpp/purify/rm_operators.h @@ -168,6 +168,27 @@ std::tuple, sopt::OperatorFunction> base_padding_an return std::make_tuple(direct, indirect); } +template +std::tuple, sopt::OperatorFunction> base_degrid_operator_1d( + const Vector &u, const Vector &weights, const t_uint imsizex, + const t_real oversample_ratio, const kernels::kernel kernel, const t_uint Ju, + const fftw_plan ft_plan) { + std::function kernelu, ftkernelu; + std::tie(ftkernelu, kernelu) = purify::create_radial_ftkernel(kernel, Ju, oversample_ratio); + sopt::OperatorFunction directFZ, indirectFZ; + std::tie(directFZ, indirectFZ) = + base_padding_and_FFT_1d(ftkernelu, imsizex, oversample_ratio, ft_plan); + sopt::OperatorFunction directG, indirectG; + PURIFY_LOW_LOG("Constructing Weighting and Gridding Operators: WG"); + PURIFY_MEDIUM_LOG("Number of channels: {}", u.size()); + std::tie(directG, indirectG) = purify::operators::init_gridding_matrix_1d( + u, weights, imsizex, oversample_ratio, kernelu, Ju); + auto direct = sopt::chained_operators(directG, directFZ); + auto indirect = sopt::chained_operators(indirectFZ, indirectG); + PURIFY_LOW_LOG("Finished consturction of RM Φ."); + return std::make_tuple(direct, indirect); +} + template std::tuple, sopt::OperatorFunction> base_degrid_operator_1d( const Vector &u, const Vector &widths, const Vector &weights, @@ -225,6 +246,23 @@ std::shared_ptr> init_degrid_operator_1d( rel_error, ft_plan); return std::make_shared>(directDegrid, M, indirectDegrid, N); } +//! Return degridding operator with no channgel averaging +template +std::shared_ptr> init_degrid_operator_1d( + const Vector &u, const Vector &weights, const t_uint imsizex, + const t_real cell, const t_real oversample_ratio = 2, + const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4, + const t_uint J_max = 30, const t_real abs_error = 1e-6, const t_real rel_error = 1e-6) { + PURIFY_LOW_LOG("Cell size in Faraday Depth is {} rad/m^2", cell); + const t_real du = rm_details::pixel_to_lambda2(cell, imsizex, oversample_ratio); + const operators::fftw_plan ft_plan = operators::fftw_plan::measure; + std::array N = {0, 1, static_cast(imsizex)}; + std::array M = {0, 1, static_cast(u.size())}; + sopt::OperatorFunction directDegrid, indirectDegrid; + std::tie(directDegrid, indirectDegrid) = purify::operators::base_degrid_operator_1d( + -u / du, weights, imsizex, oversample_ratio, kernel, Ju, ft_plan); + return std::make_shared>(directDegrid, M, indirectDegrid, N); +} } // namespace measurementoperator }; // namespace purify diff --git a/cpp/tests/rm_operator.cc b/cpp/tests/rm_operator.cc index 99520b366..aa02c13bd 100644 --- a/cpp/tests/rm_operator.cc +++ b/cpp/tests/rm_operator.cc @@ -213,3 +213,45 @@ TEST_CASE("normed operator") { } } } + +TEST_CASE("normed operator with no widths") { + const t_real oversample_ratio = 2; + const t_int power_iters = 1000; + const t_real power_tol = 1e-4; + Vector u = Vector::Random(10); + const Vector widths = Vector::Zero(u.size()); + const t_uint M = u.size(); + const Vector y = Vector::Ones(u.size()); + SECTION("kb") { + const kernels::kernel kernel = kernels::kernel::kb; + for (auto& J : {4, 5, 6, 7, 8}) { + for (auto& imsize : {128, 256, 512}) { + const t_real cell = constant::pi / imsize; + const Vector init = Vector::Ones(imsize); + auto power_method_result = sopt::algorithm::normalise_operator>( + measurementoperator::init_degrid_operator_1d>( + u * imsize / 2, Vector::Ones(M), imsize, cell, oversample_ratio, kernel, + J), + power_iters, power_tol, init); + const t_real norm = std::get<0>(power_method_result); + const std::shared_ptr>> measure_op = + std::get<2>(power_method_result); + + Vector input = Vector::Zero(imsize); + input(static_cast(imsize * 0.5)) = 1.; + const Vector y_test = + (*measure_op * input).eval() * + std::sqrt(static_cast(input.size()) * oversample_ratio); + CAPTURE(y_test.cwiseAbs().mean()); + CAPTURE(y); + CAPTURE(y_test); + CAPTURE(J); + CAPTURE(imsize) + CHECK((y_test * norm).isApprox(y, 1e-3)); + const Vector psf = + (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M * norm; + CHECK(std::real(psf(static_cast(imsize * 0.5))) == Approx(1.).margin(0.001)); + } + } + } +}