PURIFY
Next-generation radio interferometric imaging
Functions
sara_padmm_random_coverage.cc File Reference
#include "purify/config.h"
#include "purify/types.h"
#include <array>
#include <memory>
#include <random>
#include <boost/math/special_functions/erf.hpp>
#include "purify/directories.h"
#include "purify/operators.h"
#include "purify/pfitsio.h"
#include "purify/utilities.h"
#include <sopt/imaging_padmm.h>
#include <sopt/positive_quadrant.h>
#include <sopt/power_method.h>
#include <sopt/relative_variation.h>
#include <sopt/reweighted.h>
#include <sopt/utilities.h>
#include <sopt/wavelets.h>
#include <sopt/wavelets/sara.h>
+ Include dependency graph for sara_padmm_random_coverage.cc:

Go to the source code of this file.

Functions

int main (int, char **)
 

Function Documentation

◆ main()

int main ( int  ,
char **   
)

Definition at line 20 of file sara_padmm_random_coverage.cc.

20  {
21  using namespace purify;
22  sopt::logging::set_level("info");
23 
24  std::string const fitsfile = image_filename("M31.fits");
25  std::string const inputfile = output_filename("M31_input.fits");
26  std::string const outfile = output_filename("M31.tiff");
27  std::string const outfile_fits = output_filename("M31_solution.fits");
28  std::string const residual_fits = output_filename("M31_residual.fits");
29  std::string const dirty_image = output_filename("M31_dirty.tiff");
30  std::string const dirty_image_fits = output_filename("M31_dirty.fits");
31  std::string const output_vis_file = output_filename("M31_Random_coverage.vis");
32 
33  t_real const over_sample = 2;
34  auto beta = 1e-3;
35  auto M31 = pfitsio::read2d(fitsfile);
36  t_real const max = M31.array().abs().maxCoeff();
37  M31 = M31 * 1. / max;
38  pfitsio::write2d(M31.real(), inputfile);
39  // Following same formula in matlab example
40  t_real const sigma_m = constant::pi / 3;
41  // t_int const number_of_vis = std::floor(p * rho * M31.size());
42  t_int const number_of_vis = 1e4;
43  // Generating random uv(w) coverage
44  auto uv_data = utilities::random_sample_density(number_of_vis, 0, sigma_m);
45  uv_data.units = utilities::vis_units::radians;
46  utilities::write_visibility(uv_data, output_vis_file);
47  SOPT_NOTICE("Number of measurements / number of pixels: {}", uv_data.u.size() * 1. / M31.size());
48  auto measurements_transform = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
49  *measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
50  uv_data.u, uv_data.v, uv_data.w, uv_data.weights, M31.cols(), M31.rows(), over_sample,
51  kernels::kernel_from_string.at("kb"), 4, 4),
52  100, 1e-4, Vector<t_complex>::Random(M31.size())));
53 
54  sopt::wavelets::SARA const sara{
55  std::make_tuple("Dirac", 3u), std::make_tuple("DB1", 3u), std::make_tuple("DB2", 3u),
56  std::make_tuple("DB3", 3u), std::make_tuple("DB4", 3u), std::make_tuple("DB5", 3u),
57  std::make_tuple("DB6", 3u), std::make_tuple("DB7", 3u), std::make_tuple("DB8", 3u)};
58 
59  auto const Psi = sopt::linear_transform<t_complex>(sara, M31.rows(), M31.cols());
60 
61  std::mt19937_64 mersenne;
62  Vector<t_complex> const y0 =
63  (measurements_transform * Vector<t_complex>::Map(M31.data(), M31.size()));
64  // working out value of signal given SNR of 30
65  t_real sigma = utilities::SNR_to_standard_deviation(y0, 30.);
66  // adding noise to visibilities
67  uv_data.vis = utilities::add_noise(y0, 0., sigma);
68  Vector<> dimage = (measurements_transform.adjoint() * uv_data.vis).real();
69  t_real const max_val = dimage.array().abs().maxCoeff();
70  dimage = dimage / max_val;
71  Vector<t_complex> initial_estimate = Vector<t_complex>::Zero(dimage.size());
72  sopt::utilities::write_tiff(Image<t_real>::Map(dimage.data(), M31.rows(), M31.cols()),
73  dirty_image);
74  pfitsio::write2d(Image<t_real>::Map(dimage.data(), M31.rows(), M31.cols()), dirty_image_fits);
75 
76  auto const epsilon = utilities::calculate_l2_radius(uv_data.vis.size(), sigma);
77  auto const purify_regulariser_strength =
78  (Psi.adjoint() * (measurements_transform.adjoint() * uv_data.vis).eval())
79  .cwiseAbs()
80  .maxCoeff() *
81  beta;
82 
83  // auto purify_regulariser_strength = 3 * utilities::median((Psi.adjoint() *
84  // (measurements_transform.adjoint() * (uv_data.vis - y0))).real().cwiseAbs())/0.6745;
85 
86  SOPT_INFO("Using epsilon of {} \n", epsilon);
87  auto const padmm = sopt::algorithm::ImagingProximalADMM<t_complex>(uv_data.vis)
88  .itermax(1000)
89  .regulariser_strength(purify_regulariser_strength)
90  .relative_variation(1e-6)
91  .l2ball_proximal_epsilon(epsilon)
92  .tight_frame(false)
93  .l1_proximal_tolerance(1e-4)
94  .l1_proximal_nu(1)
95  .l1_proximal_itermax(50)
96  .l1_proximal_positivity_constraint(true)
97  .l1_proximal_real_constraint(true)
98  .residual_convergence(epsilon * 1.001)
99  .lagrange_update_scale(0.9)
100  .Psi(Psi)
101  .Phi(measurements_transform);
102 
103  auto const posq = sopt::algorithm::positive_quadrant(padmm);
104  auto const min_delta = sigma * std::sqrt(y0.size()) / std::sqrt(9 * M31.size());
105  // Sets weight after each padmm iteration.
106  // In practice, this means replacing the proximal of the l1 objective function.
107  auto const reweighted =
108  sopt::algorithm::reweighted(padmm).itermax(10).min_delta(min_delta).is_converged(
109  sopt::RelativeVariation<std::complex<t_real>>(1e-3));
110  auto const diagnostic = reweighted();
111  assert(diagnostic.algo.x.size() == M31.size());
112  Image<t_complex> image = Image<t_complex>::Map(diagnostic.algo.x.data(), M31.rows(), M31.cols());
113  sopt::utilities::write_tiff(image.real(), outfile);
114  pfitsio::write2d(image.real(), outfile_fits);
115  Image<t_complex> residual = Image<t_complex>::Map(
116  (measurements_transform.adjoint() *
117  (y0 - measurements_transform * Vector<t_complex>::Map(image.data(), image.size())))
118  .eval()
119  .data(),
120  M31.rows(), M31.cols());
121  pfitsio::write2d(residual.real(), residual_fits);
122  return 0;
123 }
mersenne
std::unique_ptr< std::mt19937_64 > mersenne(new std::mt19937_64(0))
purify::constant::pi
const t_real pi
mathematical constant
Definition: types.h:70
purify::kernels::kernel_from_string
const std::map< std::string, kernel > kernel_from_string
Definition: kernels.h:16
purify::logging::set_level
void set_level(const std::string &level)
Method to set the logging level of the default Log object.
Definition: logging.h:137
purify::measurementoperator::init_degrid_operator_2d
std::shared_ptr< sopt::LinearTransform< T > > init_degrid_operator_2d(const Vector< t_real > &u, const Vector< t_real > &v, const Vector< t_real > &w, const Vector< t_complex > &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 bool w_stacking=false, const t_real &cellx=1, const t_real &celly=1)
Returns linear transform that is the standard degridding operator.
Definition: operators.h:608
purify::pfitsio::write2d
void write2d(const Image< t_real > &eigen_image, const pfitsio::header_params &header, const bool &overwrite)
Write image to fits file.
Definition: pfitsio.cc:30
purify::pfitsio::read2d
Image< t_complex > read2d(const std::string &fits_name)
Read image from fits file.
Definition: pfitsio.cc:109
purify::utilities::write_visibility
void write_visibility(const utilities::vis_params &uv_vis, const std::string &file_name, const bool w_term)
Writes visibilities to txt.
Definition: uvw_utilities.cc:243
purify::utilities::SNR_to_standard_deviation
t_real SNR_to_standard_deviation(const Vector< t_complex > &y0, const t_real &SNR)
Converts SNR to RMS noise.
Definition: utilities.cc:101
purify::utilities::add_noise
Vector< t_complex > add_noise(const Vector< t_complex > &y0, const t_complex &mean, const t_real &standard_deviation)
Add guassian noise to vector.
Definition: utilities.cc:113
purify::utilities::random_sample_density
utilities::vis_params random_sample_density(const t_int vis_num, const t_real mean, const t_real standard_deviation, const t_real rms_w)
Generates a random visibility coverage.
Definition: uvw_utilities.cc:104
purify::utilities::calculate_l2_radius
t_real calculate_l2_radius(const t_uint y_size, const t_real &sigma, const t_real &n_sigma, const std::string distirbution)
A function that calculates the l2 ball radius for sopt.
Definition: utilities.cc:75
purify::output_filename
std::string output_filename(std::string const &filename)
Test output file.
Definition: directories.in.h:49
purify::image_filename
std::string image_filename(std::string const &filename)
Image filename.
Definition: directories.in.h:20
padmm
void padmm(const std::string &name, const Image< t_complex > &M31, const std::string &kernel, const t_int J, const utilities::vis_params &uv_data, const t_real sigma, const std::tuple< bool, t_real > &w_term)
Definition: padmm_random_coverage.cc:25

References purify::utilities::add_noise(), purify::utilities::calculate_l2_radius(), purify::image_filename(), purify::measurementoperator::init_degrid_operator_2d(), purify::kernels::kernel_from_string, mersenne(), purify::output_filename(), padmm(), purify::constant::pi, purify::utilities::radians, purify::utilities::random_sample_density(), purify::pfitsio::read2d(), purify::logging::set_level(), purify::utilities::SNR_to_standard_deviation(), purify::pfitsio::write2d(), and purify::utilities::write_visibility().

Generated by doxygen 1.9.1