purify::utilities Namespace Reference

Classes
struct	vis_params

Enumerations
enum class	vis_source { measurements , simulation }
enum class	vis_units { lambda , radians , pixels }

Functions
void	regroup (vis_params &uv_params, std::vector< t_int > const &groups_, const t_int max_groups)
void	regroup (vis_params &uv_params, std::vector< t_int > &image_index, std::vector< t_int > const &groups_, const t_int max_groups)
vis_params	regroup_and_scatter (vis_params const &params, std::vector< t_int > const &groups, sopt::mpi::Communicator const &comm)
std::tuple< vis_params, std::vector< t_int > >	regroup_and_all_to_all (vis_params const &params, const std::vector< t_int > &image_index, std::vector< t_int > const &groups, sopt::mpi::Communicator const &comm)
vis_params	regroup_and_all_to_all (vis_params const &params, std::vector< t_int > const &groups, sopt::mpi::Communicator const &comm)
vis_params	all_to_all_visibilities (vis_params const &params, std::vector< t_int > const &sizes, sopt::mpi::Communicator const &comm)
vis_params	scatter_visibilities (vis_params const &params, std::vector< t_int > const &sizes, sopt::mpi::Communicator const &comm)
vis_params	scatter_visibilities (sopt::mpi::Communicator const &comm)
utilities::vis_params	distribute_params (utilities::vis_params const &params, sopt::mpi::Communicator const &comm)
utilities::vis_params	set_cell_size (const sopt::mpi::Communicator &comm, utilities::vis_params const &uv_vis, const t_real &cell_x, const t_real &cell_y)
utilities::vis_params	w_stacking (utilities::vis_params const &params, sopt::mpi::Communicator const &comm, const t_int iters, const std::function< t_real(t_real)> &cost, const t_real k_means_rel_diff)
std::tuple< utilities::vis_params, std::vector< t_int >, std::vector< t_real > >	w_stacking_with_all_to_all (utilities::vis_params const &params, const t_real du, const t_int min_support, const t_int max_support, sopt::mpi::Communicator const &comm, const t_int iters, const t_real fill_relaxation, const std::function< t_real(t_real)> &cost, const t_real k_means_rel_diff)
template<class T >
t_real	step_size (T const &vis, const std::shared_ptr< sopt::LinearTransform< T > const > &measurements, const std::shared_ptr< sopt::LinearTransform< T > const > &wavelets, const t_uint sara_size)
	Calculate step size using MPI (does not include factor of 1e-3) More...
t_int	sub2ind (const t_int &row, const t_int &col, const t_int &rows, const t_int &cols)
	Converts from subscript to index for matrix. More...
std::tuple< t_int, t_int >	ind2sub (const t_int &sub, const t_int &cols, const t_int &rows)
	Converts from index to subscript for matrix. More...
t_real	mod (const t_real &x, const t_real &y)
	Mod function modified to wrap circularly for negative numbers. More...
Image< t_complex >	convolution_operator (const Image< t_complex > &a, const Image< t_complex > &b)
	Calculates the convolution between two images. More...
t_real	calculate_l2_radius (const t_uint y_size, const t_real &sigma=0, const t_real &n_sigma=2., const std::string distirbution="chi")
	A function that calculates the l2 ball radius for sopt. More...
t_real	SNR_to_standard_deviation (const Vector< t_complex > &y0, const t_real &SNR=30.)
	Converts SNR to RMS noise. More...
Vector< t_complex >	add_noise (const Vector< t_complex > &y, const t_complex &mean, const t_real &standard_deviation)
	Add guassian noise to vector. More...
bool	file_exists (const std::string &name)
	Test to see if file exists. More...
std::tuple< t_real, t_real, t_real >	fit_fwhm (const Image< t_real > &psf, const t_int &size=3)
	Method to fit Gaussian to PSF. More...
t_real	median (const Vector< t_real > &input)
	Return median of real vector. More...
t_real	dynamic_range (const Image< t_complex > &model, const Image< t_complex > &residuals, const t_real &operator_norm=1)
	Calculate the dynamic range between the model and residuals. More...
Array< t_complex >	init_weights (const Vector< t_real > &u, const Vector< t_real > &v, const Vector< t_complex > &weights, const t_real &oversample_factor, const std::string &weighting_type, const t_real &R, const t_int &ftsizeu, const t_int &ftsizev)
	Calculate weightings. More...
std::tuple< t_int, t_real >	checkpoint_log (const std::string &diagnostic)
	Reads a diagnostic file and updates parameters. More...
Matrix< t_complex >	re_sample_ft_grid (const Matrix< t_complex > &input, const t_real &re_sample_factor)
	zero pads ft grid for image up sampling and downsampling More...
Matrix< t_complex >	re_sample_image (const Matrix< t_complex > &input, const t_real &re_sample_ratio)
	resamples image size More...
template<class K >
K::Scalar	mean (const K x)
	Calculate mean of vector. More...
template<class K >
t_real	variance (const K x)
	Calculate variance of vector. More...
template<class K >
t_real	standard_deviation (const K x)
	Calculates the standard deviation of a vector. More...
template<class T0 , class T1 >
std::enable_if< std::is_same< typename T0::Scalar, typename T1::Scalar >::value and T0::IsRowMajor, Vector< typename T0::Scalar > >::type	sparse_multiply_matrix (const Eigen::SparseMatrixBase< T0 > &M, const Eigen::MatrixBase< T1 > &x)
	Parallel multiplication with a sparse matrix and vector. More...
template<class K , class L >
Image< t_complex >	parallel_multiply_image (const K &A, const L &B)
	Multiply images coefficient-wise using openmp. More...
bool	has_suffix (const std::string &str, const std::string &suff)
Matrix< t_real >	generate_antennas (const t_uint N, const t_real scale)
	Generate guassianly distributed (sigma = scale) antenna positions. More...
utilities::vis_params	antenna_to_coverage (const Matrix< t_real > &B, const t_real frequency, const t_real times, const t_real theta_ra, const t_real phi_dec, const t_real latitude)
	Provided antenna positions generate a coverage for a fixed frequency a fixed time. More...
utilities::vis_params	antenna_to_coverage (const Matrix< t_real > &B, const t_real frequency, const std::vector< t_real > &times, const t_real theta_ra, const t_real phi_dec, const t_real latitude)
	Provided antenna positions generate a coverage for a fixed frequency for mulitple times. More...
utilities::vis_params	antenna_to_coverage (const Matrix< t_real > &B, const std::vector< t_real > &frequencies, const t_real times, const t_real theta_ra, const t_real phi_dec, const t_real latitude)
	Provided antenna positions generate a coverage for multiple frequency a fixed time. More...
utilities::vis_params	antenna_to_coverage (const Matrix< t_real > &B, const std::vector< t_real > &frequencies, const std::vector< t_real > times, const t_real theta_ra, const t_real phi_dec, const t_real latitude)
Matrix< t_real >	read_ant_positions (const std::string &pos_name)
	Read in a text file of antenna positions into a matrix [x, y ,z]. More...
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=0)
	Generates a random visibility coverage. More...
utilities::vis_params	read_visibility (const std::vector< std::string > &names, const bool w_term=false)
	Read visibility files from name of vector. More...
utilities::vis_params	read_visibility (const std::string &vis_name2, const utilities::vis_params &u1)
	Reads in two visibility files. More...
t_real	streamtoreal (std::ifstream &stream)
	Reading reals from visibility file (including nan's and inf's) More...
utilities::vis_params	read_visibility_csv (const std::string &vis_name, const bool w_term=false)
	Reads in visibility csv file. More...
utilities::vis_params	read_visibility (const std::string &vis_name, const bool w_term=false)
	Reads in visibility file. More...
void	write_visibility (const utilities::vis_params &uv_vis, const std::string &file_name, const bool w_term=false)
	Writes visibilities to txt. More...
utilities::vis_params	set_cell_size (const utilities::vis_params &uv_vis, const t_real &max_u, const t_real &max_v, const t_real &input_cell_size_u, const t_real &input_cell_size_v)
utilities::vis_params	set_cell_size (const utilities::vis_params &uv_vis, const t_real &cell_size_u=0, const t_real &cell_size_v=0)
	Scales visibilities to a given pixel size in arcseconds. More...
utilities::vis_params	uv_scale (const utilities::vis_params &uv_vis, const t_int &ftsizeu, const t_int &ftsizev)
	scales the visibilities to units of pixels More...
utilities::vis_params	convert_to_pixels (const utilities::vis_params &uv_vis, const t_real cell_x, const t_real cell_y, const t_real imsizex, const t_real imsizey, const t_real oversample_ratio)
	Converts u and v coordaintes to units of pixels. More...
utilities::vis_params	conjugate_w (const utilities::vis_params &uv_vis)
	reflects visibilities into the w >= 0 domain More...
template<class T >
T	calculate_rotated_u (const T &u, const T &v, const T &w, const t_real alpha, const t_real beta, const t_real gamma)
	calculate the rotated u from euler angles in zyz and starting coordinates (u, v, w) More...
template<class T >
T	calculate_rotated_v (const T &u, const T &v, const T &w, const t_real alpha, const t_real beta, const t_real gamma)
	calculate the rotated v from euler angles in zyz and starting coordinates (u, v, w) More...
template<class T >
T	calculate_rotated_w (const T &u, const T &v, const T &w, const t_real alpha, const t_real beta, const t_real gamma)
	calculate the rotated w from euler angles in zyz and starting coordinates (u, v, w) More...
template<class T >
utilities::vis_params	antenna_to_coverage (const t_uint N, const t_real scale, const T &frequency)
utilities::vis_params	antenna_to_coverage (const t_uint N)
	Using guassianly distributed (sigma = pi) antenna positions generate a coverage. More...
template<class F >
utilities::vis_params	antenna_to_coverage_general (const Matrix< t_real > &B, const std::vector< t_real > &frequencies, const std::vector< t_real > &times, const F &position_angle_RA_function, const F &position_angle_DEC_function, const t_real latitude)
template<class T , class K >
utilities::vis_params	antenna_to_coverage (const Vector< t_real > &x, const Vector< t_real > &y, const Vector< t_real > &z, const T &frequencies, const K &times, const t_real theta_ra, const t_real phi_dec, const t_real latitude)
	Provided antenna positions generate a coverage. More...
template<class T , class K >
utilities::vis_params	read_ant_positions_to_coverage (const std::string &pos_name, const T &frequencies, const K &times, const t_real theta_ra, const t_real phi_dec, const t_real latitude)
vis_params	sort_by_w (const vis_params &uv_data)
	Sort visibilities to be from w_max to w_min. More...
template<typename T >
t_real	w_lambert (T x, const t_real &relative_difference)
	Calculate W Lambert function. More...
template<typename T >
t_real	w_lambert (const T &x, const t_real &relative_difference=1e-8)

Enumeration Type Documentation

vis_source

enum purify::utilities::vis_source

strong

Enumerator
measurements
simulation

Definition at line 13 of file uvw_utilities.h.

{ measurements, simulation };

vis_units

enum purify::utilities::vis_units

strong

Enumerator
lambda
radians
pixels

Definition at line 14 of file uvw_utilities.h.

{ lambda, radians, pixels };

Function Documentation

add_noise()

Vector< t_complex > purify::utilities::add_noise	(	const Vector< t_complex > &	y0,
		const t_complex &	mean,
		const t_real &	standard_deviation
	)

Add guassian noise to vector.

Definition at line 113 of file utilities.cc.

                                                              {
  /*
          Adds complex valued Gaussian noise to vector
          y0:: vector before noise
          mean:: complex valued mean of noise
          standard_deviation:: rms of noise
  */
  auto sample = [&mean, &standard_deviation](t_complex x) {
    std::random_device rd_real;
    std::random_device rd_imag;
    std::mt19937_64 rng_real(rd_real());
    std::mt19937_64 rng_imag(rd_imag());
    static std::normal_distribution<t_real> normal_dist_real(std::real(mean), standard_deviation);
    static std::normal_distribution<t_real> normal_dist_imag(std::imag(mean), standard_deviation);
    t_complex I(0, 1.);
    return normal_dist_real(rng_real) + I * normal_dist_imag(rng_imag);
  };
 
  auto noise = Vector<t_complex>::Zero(y0.size()).unaryExpr(sample);
 
  return y0 + noise;
}

References purify::I, mean(), and standard_deviation().

Referenced by b_utilities::dirty_measurements(), dirty_visibilities(), getInputData(), and main().

all_to_all_visibilities()

vis_params purify::utilities::all_to_all_visibilities	(	vis_params const &	params,
		std::vector< t_int > const &	sizes,
		sopt::mpi::Communicator const &	comm
	)

Definition at line 122 of file mpi_utilities.cc.

                                                                    {
  if (comm.size() == 1) return params;
  vis_params result;
  result.u = comm.all_to_allv(params.u, sizes);
  result.v = comm.all_to_allv(params.v, sizes);
  result.w = comm.all_to_allv(params.w, sizes);
  result.vis = comm.all_to_allv(params.vis, sizes);
  result.weights = comm.all_to_allv(params.weights, sizes);
  result.units = static_cast<utilities::vis_units>(comm.broadcast(static_cast<int>(params.units)));
  result.ra = comm.broadcast(params.ra);
  result.dec = comm.broadcast(params.dec);
  result.average_frequency = comm.broadcast(params.average_frequency);
  return result;
}

References purify::utilities::vis_params::average_frequency, purify::utilities::vis_params::dec, purify::utilities::vis_params::ra, purify::utilities::vis_params::u, purify::utilities::vis_params::units, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

Referenced by regroup_and_all_to_all().

antenna_to_coverage() [1/7]

utilities::vis_params purify::utilities::antenna_to_coverage	(	const Matrix< t_real > &	B,
		const std::vector< t_real > &	frequencies,
		const std::vector< t_real >	times,
		const t_real	theta_ra,
		const t_real	phi_dec,
		const t_real	latitude
	)

Provided antenna positions generate a coverage for multiple frequencies for multiple times for earth rotation synthesis

Definition at line 65 of file uvw_utilities.cc.

                                                                                       {
  auto const position_angle_RA_function = [theta_ra](const t_real t) -> t_real {
    return theta_ra + constant::omega_e * t;
  };
  auto const position_angle_DEC_function = [phi_dec](const t_real t) -> t_real { return phi_dec; };
  return antenna_to_coverage_general<std::function<t_real(t_real)>>(
      B, frequencies, times, position_angle_RA_function, position_angle_DEC_function, latitude);
}

References antenna_to_coverage_general(), and purify::constant::omega_e.

antenna_to_coverage() [2/7]

utilities::vis_params purify::utilities::antenna_to_coverage	(	const Matrix< t_real > &	B,
		const std::vector< t_real > &	frequency,
		const t_real	times,
		const t_real	theta_ra,
		const t_real	phi_dec,
		const t_real	latitude
	)

Provided antenna positions generate a coverage for multiple frequency a fixed time.

Definition at line 57 of file uvw_utilities.cc.

                                                                 {
  return antenna_to_coverage(B, frequency, std::vector<t_real>(1, times), theta_ra, phi_dec,
                             latitude);
}

References antenna_to_coverage().

antenna_to_coverage() [3/7]

utilities::vis_params purify::utilities::antenna_to_coverage	(	const Matrix< t_real > &	B,
		const t_real	frequency,
		const std::vector< t_real > &	times,
		const t_real	theta_ra,
		const t_real	phi_dec,
		const t_real	latitude
	)

Provided antenna positions generate a coverage for a fixed frequency for mulitple times.

Definition at line 50 of file uvw_utilities.cc.

                                                                                       {
  return antenna_to_coverage(B, std::vector<t_real>(1, frequency), times, theta_ra, phi_dec,
                             latitude);
}

References antenna_to_coverage().

antenna_to_coverage() [4/7]

utilities::vis_params purify::utilities::antenna_to_coverage	(	const Matrix< t_real > &	B,
		const t_real	frequency,
		const t_real	times,
		const t_real	theta_ra,
		const t_real	phi_dec,
		const t_real	latitude
	)

Provided antenna positions generate a coverage for a fixed frequency a fixed time.

Definition at line 43 of file uvw_utilities.cc.

                                                                                       {
  return antenna_to_coverage(B, std::vector<t_real>(1, frequency), std::vector<t_real>(1, times),
                             theta_ra, phi_dec, latitude);
}

Referenced by antenna_to_coverage(), read_ant_positions_to_coverage(), and TEST_CASE().

antenna_to_coverage() [5/7]

utilities::vis_params purify::utilities::antenna_to_coverage

(

const t_uint

N

)

Using guassianly distributed (sigma = pi) antenna positions generate a coverage.

antenna_to_coverage() [6/7]

template<class T >

utilities::vis_params purify::utilities::antenna_to_coverage	(	const t_uint	N,
		const t_real	scale,
		const T &	frequency
	)

Definition at line 92 of file uvw_utilities.h.

                                                                                                {
  return antenna_to_coverage(generate_antennas(N, scale), frequency, 0., 0., 0., 0.);
}

References antenna_to_coverage(), and generate_antennas().

antenna_to_coverage() [7/7]

template<class T , class K >

utilities::vis_params purify::utilities::antenna_to_coverage	(	const Vector< t_real > &	x,
		const Vector< t_real > &	y,
		const Vector< t_real > &	z,
		const T &	frequencies,
		const K &	times,
		const t_real	theta_ra,
		const t_real	phi_dec,
		const t_real	latitude
	)

Provided antenna positions generate a coverage.

Definition at line 166 of file uvw_utilities.h.

                                                                                       {
  if (x.size() != y.size()) throw std::runtime_error("x and y positions are not the same amount.");
  if (x.size() != z.size()) throw std::runtime_error("x and z positions are not the same amount.");
  Matrix<t_real> B(x.size(), 3);
  B.col(0) = x;
  B.col(1) = y;
  B.col(2) = z;
  return antenna_to_coverage(B, frequencies, times, theta_ra, phi_dec, latitude);
}

References antenna_to_coverage().

antenna_to_coverage_general()

template<class F >

utilities::vis_params purify::utilities::antenna_to_coverage_general	(	const Matrix< t_real > &	B,
		const std::vector< t_real > &	frequencies,
		const std::vector< t_real > &	times,
		const F &	position_angle_RA_function,
		const F &	position_angle_DEC_function,
		const t_real	latitude
	)

Provided antenna positions generate a coverage for multiple frequencies for multiple times at different pointings

Definition at line 119 of file uvw_utilities.h.

                                                                         {
  if (B.cols() != 3) throw std::runtime_error("Antennae coordinates are not 3D vectors.");
  const t_uint M = B.rows() * (B.rows() - 1) / 2;
  const t_uint chans = frequencies.size();
  const t_uint time_samples = times.size();
  Vector<t_real> u = Vector<t_real>::Zero(M * chans * time_samples);
  Vector<t_real> v = Vector<t_real>::Zero(M * chans * time_samples);
  Vector<t_real> w = Vector<t_real>::Zero(M * chans * time_samples);
  Vector<t_complex> weights = Vector<t_complex>::Ones(M * chans * time_samples);
  Vector<t_complex> vis = Vector<t_complex>::Ones(M * chans * time_samples);
  t_uint m = 0;
  for (t_uint i = 0; i < B.rows(); i++) {
    for (t_uint j = i + 1; j < B.rows(); j++) {
      const Vector<t_real> r = (B.row(i) - B.row(j));
      for (t_int k = 0; k < chans; k++) {
        for (t_int t = 0; t < time_samples; t++) {
          const t_int index = m + M * (k + chans * t);
          u(index) = utilities::calculate_rotated_u(
                         r(0), r(1), r(2), position_angle_RA_function(times.at(t)),
                         position_angle_DEC_function(times.at(t)) - latitude, 0.) *
                     frequencies.at(k) / constant::c;
          v(index) = utilities::calculate_rotated_v(
                         r(0), r(1), r(2), position_angle_RA_function(times.at(t)),
                         position_angle_DEC_function(times.at(t)) - latitude, 0.) *
                     frequencies.at(k) / constant::c;
          w(index) = utilities::calculate_rotated_w(
                         r(0), r(1), r(2), position_angle_RA_function(times.at(t)),
                         position_angle_DEC_function(times.at(t)) - latitude, 0.) *
                     frequencies.at(k) / constant::c;
        }
      }
      m++;
    }
  }
  if (M != m)
    throw std::runtime_error(
        "Number of created baselines does not match expected baseline number N * (N - 1) / 2.");
  utilities::vis_params coverage(u, v, w, vis, weights, utilities::vis_units::lambda);
  return coverage;
};

References purify::constant::c, calculate_rotated_u(), calculate_rotated_v(), calculate_rotated_w(), lambda, operators_test::u, and operators_test::v.

Referenced by antenna_to_coverage().

calculate_l2_radius()

t_real purify::utilities::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 at line 75 of file utilities.cc.

                                                         {
  /*
                  Calculates the epsilon, the radius of the l2_ball in sopt
                  y:: vector for the l2 ball
  */
  if (distirbution == "chi^2") {
    if (sigma == 0) {
      return std::sqrt(2 * y_size + n_sigma * std::sqrt(4 * y_size));
    }
    return std::sqrt(2 * y_size + n_sigma * std::sqrt(4 * y_size)) * sigma;
  }
  if (distirbution == "chi") {
    auto alpha =
        1. / (8 * y_size) - 1. / (128 * y_size * y_size);  // series expansion for gamma relation
    auto mean = std::sqrt(2) * std::sqrt(y_size) *
                (1 - alpha);  // using Gamma(k+1/2)/Gamma(k) asymptotic formula
    auto standard_deviation = std::sqrt(2 * y_size * alpha * (2 - alpha));
    if (sigma == 0) {
      return mean + n_sigma * standard_deviation;
    }
    return (mean + n_sigma * standard_deviation) * sigma;
  }
 
  return 1.;
}

References mean(), and standard_deviation().

Referenced by main(), padmm(), padmm_factory(), TEST_CASE(), and b_utilities::updateMeasurements().

calculate_rotated_u()

template<class T >

T purify::utilities::calculate_rotated_u	(	const T &	u,
		const T &	v,
		const T &	w,
		const t_real	alpha,
		const t_real	beta,
		const t_real	gamma
	)

calculate the rotated u from euler angles in zyz and starting coordinates (u, v, w)

Definition at line 59 of file uvw_utilities.h.

                                          {
  return u * (std::cos(alpha) * std::cos(beta) * std::cos(gamma) -
              std::sin(alpha) * std::sin(gamma)) +
         v * (-std::cos(alpha) * std::cos(beta) * std::sin(gamma) -
              std::sin(alpha) * std::cos(gamma)) +
         w * std::cos(alpha) * std::sin(beta);
}

References operators_test::u, and operators_test::v.

Referenced by antenna_to_coverage_general(), and TEST_CASE().

calculate_rotated_v()

template<class T >

T purify::utilities::calculate_rotated_v	(	const T &	u,
		const T &	v,
		const T &	w,
		const t_real	alpha,
		const t_real	beta,
		const t_real	gamma
	)

calculate the rotated v from euler angles in zyz and starting coordinates (u, v, w)

Definition at line 69 of file uvw_utilities.h.

                                          {
  return u * (std::sin(alpha) * std::cos(beta) * std::cos(gamma) +
              std::cos(alpha) * std::sin(gamma)) +
         v * (-std::sin(alpha) * std::cos(beta) * std::sin(gamma) +
              std::cos(alpha) * std::cos(gamma)) +
         w * std::sin(alpha) * std::sin(beta);
}

References operators_test::u, and operators_test::v.

Referenced by antenna_to_coverage_general(), and TEST_CASE().

calculate_rotated_w()

template<class T >

T purify::utilities::calculate_rotated_w	(	const T &	u,
		const T &	v,
		const T &	w,
		const t_real	alpha,
		const t_real	beta,
		const t_real	gamma
	)

calculate the rotated w from euler angles in zyz and starting coordinates (u, v, w)

Definition at line 79 of file uvw_utilities.h.

                                          {
  return u * (-std::sin(beta) * std::cos(gamma)) + v * (std::sin(beta) * std::sin(gamma)) +
         w * std::cos(beta);
}

References operators_test::u, and operators_test::v.

Referenced by antenna_to_coverage_general(), and TEST_CASE().

checkpoint_log()

std::tuple< t_int, t_real > purify::utilities::checkpoint_log

(

const std::string &

diagnostic

)

Reads a diagnostic file and updates parameters.

Definition at line 287 of file utilities.cc.

                                                                    {
  // reads a log file and returns the latest parameters
  if (!utilities::file_exists(diagnostic)) return std::make_tuple(0, 0);
 
  t_int iters = 0;
  t_real gamma = 0;
  std::ifstream log_file(diagnostic);
  std::string entry;
  std::string s;
 
  std::getline(log_file, s);
 
  while (log_file) {
    if (!std::getline(log_file, s)) break;
    std::istringstream ss(s);
    std::getline(ss, entry, ' ');
    iters = std::floor(std::stod(entry));
    std::getline(ss, entry, ' ');
    gamma = std::stod(entry);
  }
  log_file.close();
  return std::make_tuple(iters, gamma);
}

References file_exists().

conjugate_w()

utilities::vis_params purify::utilities::conjugate_w

(

const utilities::vis_params &

uv_vis

)

reflects visibilities into the w >= 0 domain

Definition at line 352 of file uvw_utilities.cc.

                                                                   {
  utilities::vis_params output = uv_vis;
#pragma omp parallel for
  for (t_uint i = 0; i < output.size(); i++) {
    if (uv_vis.w(i) < 0) {
      output.w(i) = -uv_vis.w(i);
      output.v(i) = -uv_vis.v(i);
      output.u(i) = -uv_vis.u(i);
      output.vis(i) = std::conj(uv_vis.vis(i));
    }
    assert(output.w(i) >= 0);
  }
  return output;
}

References purify::utilities::vis_params::size(), purify::utilities::vis_params::u, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, and purify::utilities::vis_params::w.

Referenced by getInputData(), main(), and TEST_CASE().

convert_to_pixels()

utilities::vis_params purify::utilities::convert_to_pixels	(	const utilities::vis_params &	uv_vis,
		const t_real	cell_x,
		const t_real	cell_y,
		const t_real	imsizex,
		const t_real	imsizey,
		const t_real	oversample_ratio
	)

Converts u and v coordaintes to units of pixels.

Definition at line 332 of file uvw_utilities.cc.

                                                                                             {
  t_real du = 1;
  t_real dv = 1;
  if (uv_vis.units == utilities::vis_units::lambda) {
    du = widefield::pixel_to_lambda(cell_x, imsizex, oversample_ratio);
    dv = widefield::pixel_to_lambda(cell_y, imsizey, oversample_ratio);
  }
  if (uv_vis.units == utilities::vis_units::radians) {
    du = constant::pi / std::floor(imsizex * oversample_ratio);
    dv = constant::pi / std::floor(imsizey * oversample_ratio);
  }
  auto out = uv_vis;
  out.u = uv_vis.u / du;
  out.v = uv_vis.v / dv;
  out.units = utilities::vis_units::pixels;
  return out;
}

References lambda, purify::constant::pi, purify::widefield::pixel_to_lambda(), pixels, radians, purify::utilities::vis_params::u, purify::utilities::vis_params::units, and purify::utilities::vis_params::v.

Referenced by purify::measurementoperator::init_degrid_operator_2d().

convolution_operator()

Image< t_complex > purify::utilities::convolution_operator	(	const Image< t_complex > &	a,
		const Image< t_complex > &	b
	)

Calculates the convolution between two images.

Definition at line 50 of file utilities.cc.

                                                                                            {
  /*
  returns the convolution of images a with images b
  a:: vector a, which is shifted
  b:: vector b, which is fixed
  */
 
  // size of a image
  t_int a_y = a.rows();
  t_int a_x = a.cols();
  // size of b image
  t_int b_y = b.rows();
  t_int b_x = b.cols();
 
  Image<t_complex> output = Image<t_complex>::Zero(a_y + b_y, a_x + b_x);
 
  for (t_int l = 0; l < b.cols(); ++l) {
    for (t_int k = 0; k < b.rows(); ++k) {
      output(k, l) = (a * b.block(k, l, a_y, a_x)).sum();
    }
  }
 
  return output;
}

distribute_params()

utilities::vis_params purify::utilities::distribute_params	(	utilities::vis_params const &	params,
		sopt::mpi::Communicator const &	comm
	)

Definition at line 176 of file mpi_utilities.cc.

                                                                         {
  if (comm.is_root() and comm.size() > 1) {
    auto const order = distribute::distribute_measurements(params, comm, distribute::plan::radial);
    return utilities::regroup_and_scatter(params, order, comm);
  } else if (comm.size() > 1)
    return utilities::scatter_visibilities(comm);
  return params;
}

References purify::distribute::distribute_measurements(), purify::distribute::radial, regroup_and_scatter(), and scatter_visibilities().

dynamic_range()

t_real purify::utilities::dynamic_range	(	const Image< t_complex > &	model,
		const Image< t_complex > &	residuals,
		const t_real &	operator_norm
	)

Calculate the dynamic range between the model and residuals.

Definition at line 232 of file utilities.cc.

                                                  {
  /*
  Returns value of noise rms given a measurement vector and signal to noise ratio
  y0:: complex valued vector before noise added
  SNR:: signal to noise ratio
 
  This calculation follows Carrillo et al. (2014), PURIFY a new approach to radio interferometric
  imaging
  */
  return std::sqrt(model.size()) * (operator_norm * operator_norm) / residuals.matrix().norm() *
         model.cwiseAbs().maxCoeff();
}

file_exists()

bool purify::utilities::file_exists

(

const std::string &

name

)

Test to see if file exists.

Definition at line 137 of file utilities.cc.

                                        {
  /*
          Checks if a file exists
          name:: path of file checked for existance.
 
          returns true if exists and false if not exists
  */
  struct stat buffer;
  return (stat(name.c_str(), &buffer) == 0);
}

Referenced by checkpoint_log(), and TEST_CASE().

fit_fwhm()

std::tuple< t_real, t_real, t_real > purify::utilities::fit_fwhm	(	const Image< t_real > &	psf,
		const t_int &	size
	)

Method to fit Gaussian to PSF.

Definition at line 148 of file utilities.cc.

                                                                                       {
  /*
          Find FWHM of point spread function, using least squares.
 
          psf:: point spread function, assumed to be normalised.
 
          The method assumes that you have sampled at least
          3 pixels across the beam.
  */
 
  t_int x0 = std::floor(psf.cols() * 0.5);
  t_int y0 = std::floor(psf.rows() * 0.5);
 
  // finding patch
  Image<t_real> patch =
      psf.block(std::floor(y0 - size * 0.5) + 1, std::floor(x0 - size * 0.5) + 1, size, size);
  PURIFY_LOW_LOG("Fitting to Patch");
 
  // finding values for least squares
 
  std::vector<t_tripletList> entries;
  t_int total_entries = 0;
 
  for (t_int i = 0; i < patch.cols(); ++i) {
    for (t_int j = 0; j < patch.rows(); ++j) {
      entries.emplace_back(j, i, std::log(patch(j, i)));
      total_entries++;
    }
  }
  Matrix<t_real> A = Matrix<t_real>::Zero(total_entries, 4);
  Vector<t_real> q = Vector<t_real>::Zero(total_entries);
  // putting values into a vector and matrix for least squares
  for (t_int i = 0; i < total_entries; ++i) {
    t_real x = entries.at(i).col() - size * 0.5 + 0.5;
    t_real y = entries.at(i).row() - size * 0.5 + 0.5;
    A(i, 0) = static_cast<t_real>(x * x);  // x^2
    A(i, 1) = static_cast<t_real>(x * y);  // x * y
    A(i, 2) = static_cast<t_real>(y * y);  // y^2
    q(i) = std::real(entries.at(i).value());
  }
  // least squares fitting
  const auto solution =
      static_cast<Vector<t_real>>(A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(q));
  t_real const a = -solution(0);
  t_real const b = +solution(1) * 0.5;
  t_real const c = -solution(2);
  // parameters of Gaussian
  t_real theta = std::atan2(b, std::sqrt(std::pow(2 * b, 2) + std::pow(c - a, 2)) +
                                   (c - a) * 0.5);  // some relatively complicated trig identity to
                                                    // go from tan(2theta) to tan(theta).
  t_real t = 0.;
  t_real s = 0.;
  if (std::abs(b) < 1e-13) {
    t = a;
    s = c;
    theta = 0;
  } else {
    t = (a + c - 2 * b / std::sin(2 * theta)) * 0.5;
    s = (a + c + 2 * b / std::sin(2 * theta)) * 0.5;
  }
 
  t_real const sigma_x = std::sqrt(1 / (2 * t));
  t_real const sigma_y = std::sqrt(1 / (2 * s));
 
  std::tuple<t_real, t_real, t_real> fit_parameters;
 
  // fit for the beam rms, used for FWHM
  std::get<0>(fit_parameters) = sigma_x;
  std::get<1>(fit_parameters) = sigma_y;
  std::get<2>(fit_parameters) = theta;  // because 2*theta is periodic with pi.
  return fit_parameters;
}

References purify::constant::c, and PURIFY_LOW_LOG.

Referenced by TEST_CASE().

generate_antennas()

Matrix< t_real > purify::utilities::generate_antennas	(	const t_uint	N,
		const t_real	scale
	)

Generate guassianly distributed (sigma = scale) antenna positions.

Definition at line 19 of file uvw_utilities.cc.

                                                                     {
  Matrix<t_real> B = Matrix<t_real>::Zero(N, 3);
  const t_real mean = 0;
  const t_real standard_deviation = scale;
  auto sample = [&mean, &standard_deviation]() {
    std::random_device rd;
    std::mt19937_64 rng(rd());
    t_real output = 4 * standard_deviation + mean;
    static std::normal_distribution<> normal_dist(mean, standard_deviation);
    // ensures that all sample points are within bounds
    while (std::abs(output - mean) > 3 * standard_deviation) {
      output = normal_dist(rng);
    }
    if (output != output) PURIFY_DEBUG("New random-sample density: {}", output);
    return output;
  };
  for (t_uint i = 0; i < N; i++) {
    B(i, 0) = sample();
    B(i, 1) = sample();
    B(i, 2) = sample();
  }
  return B * scale;
}

References mean(), PURIFY_DEBUG, and standard_deviation().

Referenced by antenna_to_coverage(), and TEST_CASE().

has_suffix()

bool purify::utilities::has_suffix	(	const std::string &	str,
		const std::string &	suff
	)

Definition at line 15 of file uvw_utilities.cc.

                                                             {
  return str.size() >= suff.size() && str.compare(str.size() - suff.size(), suff.size(), suff) == 0;
}

Referenced by read_visibility().

ind2sub()

std::tuple< t_int, t_int > purify::utilities::ind2sub	(	const t_int &	sub,
		const t_int &	cols,
		const t_int &	rows
	)

Converts from index to subscript for matrix.

Definition at line 26 of file utilities.cc.

                                                                                       {
  /*
    Converts index of a matrix to (row, column). This does the same as the matlab funciton sub2ind,
    converts subscript to index.
 
    sub:: subscript of entry in matrix
    cols:: number of columns for matrix
    rows:: number of rows for matrix
    row:: output row of matrix
    col:: output column of matrix
 
   */
  return std::make_tuple<t_int, t_int>(std::floor((sub - (sub % cols)) / cols), sub % cols);
}

Referenced by purify::wproj_utilities::row_wise_convolution(), and purify::wproj_utilities::row_wise_sparse_convolution().

init_weights()

Array< t_complex > purify::utilities::init_weights	(	const Vector< t_real > &	u,
		const Vector< t_real > &	v,
		const Vector< t_complex > &	weights,
		const t_real &	oversample_factor,
		const std::string &	weighting_type,
		const t_real &	R,
		const t_int &	ftsizeu,
		const t_int &	ftsizev
	)

Calculate weightings.

Definition at line 246 of file utilities.cc.

                                                                          {
  /*
    Calculate the weights to be applied to the visibilities in the measurement operator.
    It does none, whiten, natural, uniform, and robust.
  */
  if (weighting_type == "none") return weights.array() * 0 + 1.;
  if (weighting_type == "natural" or weighting_type == "whiten") return weights;
  Vector<t_complex> out_weights(weights.size());
  if ((weighting_type == "uniform") or (weighting_type == "robust")) {
    t_real scale =
        1. / oversample_factor;  // scale for fov, controlling the region of sidelobe supression
    Matrix<t_complex> gridded_weights = Matrix<t_complex>::Zero(ftsizev, ftsizeu);
    for (t_int i = 0; i < weights.size(); ++i) {
      t_int q = utilities::mod(floor(u(i) * scale), ftsizeu);
      t_int p = utilities::mod(floor(v(i) * scale), ftsizev);
      gridded_weights(p, q) += 1 * 1;  // I get better results assuming all the weights are the
                                       // same. It looks like miriad does this as well.
    }
    t_complex const sum_weights = (weights.array() * weights.array()).sum();
    t_complex const sum_grid_weights2 = (gridded_weights.array() * gridded_weights.array()).sum();
    t_complex const robust_scale =
        sum_weights / sum_grid_weights2 * 25. *
        std::pow(10, -2 * R);  // Following standard formula, a bit different from miriad.
    gridded_weights = gridded_weights.array().sqrt();
    for (t_int i = 0; i < weights.size(); ++i) {
      t_int q = utilities::mod(floor(u(i) * scale), ftsizeu);
      t_int p = utilities::mod(floor(v(i) * scale), ftsizev);
      if (weighting_type == "uniform") out_weights(i) = weights(i) / gridded_weights(p, q);
      if (weighting_type == "robust") {
        out_weights(i) = weights(i) / std::sqrt(1. + robust_scale * gridded_weights(p, q) *
                                                         gridded_weights(p, q));
      }
    }
  } else
    throw std::runtime_error("Wrong weighting type, " + weighting_type + " not recognised.");
  return out_weights;
}

References mod(), operators_test::u, and operators_test::v.

mean()

template<class K >

K::Scalar purify::utilities::mean

(

const K

x

)

Calculate mean of vector.

Definition at line 21 of file utilities.h.

                                 {
  // Calculate mean of vector x
  return x.array().mean();
}

Referenced by add_noise(), calculate_l2_radius(), generate_antennas(), random_sample_density(), and TEST_CASE().

median()

t_real purify::utilities::median

(

const Vector< t_real > &

input

)

Return median of real vector.

Definition at line 221 of file utilities.cc.

                                           {
  // Finds the median of a real vector x
  auto size = input.size();
  Vector<t_real> x(size);
  std::copy(input.data(), input.data() + size, x.data());
  std::sort(x.data(), x.data() + size);
  if (std::floor(size / 2) - std::ceil(size / 2) == 0)
    return 0.5 * (x(int(std::floor(size / 2) - 1)) + x(int(std::floor(size / 2))));
  return x(int(std::ceil(size / 2)));
}

Referenced by TEST_CASE().

mod()

t_real purify::utilities::mod	(	const t_real &	x,
		const t_real &	y
	)

Mod function modified to wrap circularly for negative numbers.

Definition at line 41 of file utilities.cc.

                                             {
  /*
    returns x % y, and warps circularly around y for negative values.
  */
  t_real r = std::fmod(x, y);
  if (r < 0) r = y + r;
  return r;
}

Referenced by purify::widefield::estimate_sample_density(), purify::details::init_gridding_matrix_2d(), purify::operators::init_on_the_fly_gridding_matrix_2d(), init_weights(), purify::wproj_utilities::row_wise_convolution(), purify::wproj_utilities::row_wise_sparse_convolution(), purify::widefield::sample_density_weights(), and TEST_CASE().

parallel_multiply_image()

template<class K , class L >

Image<t_complex> purify::utilities::parallel_multiply_image	(	const K &	A,
		const L &	B
	)

Multiply images coefficient-wise using openmp.

Definition at line 83 of file utilities.h.

                                                                 {
  const t_int rows = A.rows();
  const t_int cols = A.cols();
  Image<t_complex> C = Matrix<t_complex>::Zero(rows, cols);
 
#pragma omp simd collapse(2)
  for (t_int i = 0; i < cols; ++i)
    for (t_int j = 0; j < rows; ++j) C(j, i) = A(j, i) * B(j, i);
  return C;
}

random_sample_density()

utilities::vis_params purify::utilities::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 at line 104 of file uvw_utilities.cc.

                                                                                                 {
  /*
          Generates a random sampling density for visibility coverage
          vis_num:: number of visibilities
          mean:: mean of distribution
          standard_deviation:: standard deviation of distirbution
  */
  auto sample = [&mean, &standard_deviation]() {
    std::random_device rd;
    std::mt19937_64 rng(rd());
    t_real output = 4 * standard_deviation + mean;
    static std::normal_distribution<> normal_dist(mean, standard_deviation);
    // ensures that all sample points are within bounds
    while (std::abs(output - mean) > 3 * standard_deviation) {
      output = normal_dist(rng);
    }
    if (output != output) PURIFY_DEBUG("New random-sample density: {}", output);
    return output;
  };
 
  utilities::vis_params uv_vis;
  uv_vis.u = Vector<t_real>::Zero(vis_num);
  uv_vis.v = Vector<t_real>::Zero(vis_num);
  uv_vis.w = Vector<t_real>::Zero(vis_num);
  // #pragma omp parallel for
  for (t_uint i = 0; i < vis_num; i++) {
    uv_vis.u(i) = sample();
    uv_vis.v(i) = sample();
    uv_vis.w(i) = sample();
  }
  uv_vis.w = uv_vis.w / standard_deviation * rms_w;
  uv_vis.weights = Vector<t_complex>::Constant(vis_num, 1);
  uv_vis.vis = Vector<t_complex>::Constant(vis_num, 1);
  uv_vis.ra = 0;
  uv_vis.dec = 0;
  uv_vis.average_frequency = 0;
  return uv_vis;
}

References purify::utilities::vis_params::average_frequency, purify::utilities::vis_params::dec, mean(), PURIFY_DEBUG, purify::utilities::vis_params::ra, standard_deviation(), purify::utilities::vis_params::u, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

Referenced by dirty_visibilities(), getInputData(), main(), b_utilities::random_measurements(), GridOperatorFixture::SetUp(), DegridOperatorFixture::SetUp(), and TEST_CASE().

re_sample_ft_grid()

Matrix< t_complex > purify::utilities::re_sample_ft_grid	(	const Matrix< t_complex > &	input,
		const t_real &	re_sample_ratio
	)

zero pads ft grid for image up sampling and downsampling

Definition at line 311 of file utilities.cc.

                                                                                                   {
  /*
    up samples image using by zero padding the fft
 
    input:: fft to be upsampled, with zero frequency at (0,0) of the matrix.
 
  */
  if (re_sample_ratio == 1) return input;
 
  // sets up dimensions for old image
  t_int old_x_centre_floor = std::floor(input.cols() * 0.5);
  t_int old_y_centre_floor = std::floor(input.rows() * 0.5);
  // need ceilling in case image is of odd dimension
  t_int old_x_centre_ceil = std::ceil(input.cols() * 0.5);
  t_int old_y_centre_ceil = std::ceil(input.rows() * 0.5);
 
  // sets up dimensions for new image
  t_int new_x = std::floor(input.cols() * re_sample_ratio);
  t_int new_y = std::floor(input.rows() * re_sample_ratio);
 
  t_int new_x_centre_floor = std::floor(new_x * 0.5);
  t_int new_y_centre_floor = std::floor(new_y * 0.5);
  // need ceilling in case image is of odd dimension
  t_int new_x_centre_ceil = std::ceil(new_x * 0.5);
  t_int new_y_centre_ceil = std::ceil(new_y * 0.5);
 
  Matrix<t_complex> output = Matrix<t_complex>::Zero(new_y, new_x);
 
  t_int box_dim_x;
  t_int box_dim_y;
 
  // now have to move each quadrant into new grid
  box_dim_x = std::min(old_x_centre_ceil, new_x_centre_ceil);
  box_dim_y = std::min(old_y_centre_ceil, new_y_centre_ceil);
  //(0, 0)
  output.bottomLeftCorner(box_dim_y, box_dim_x) = input.bottomLeftCorner(box_dim_y, box_dim_x);
 
  box_dim_x = std::min(old_x_centre_ceil, new_x_centre_ceil);
  box_dim_y = std::min(old_y_centre_floor, new_y_centre_floor);
  //(y0, 0)
  output.topLeftCorner(box_dim_y, box_dim_x) = input.topLeftCorner(box_dim_y, box_dim_x);
 
  box_dim_x = std::min(old_x_centre_floor, new_x_centre_floor);
  box_dim_y = std::min(old_y_centre_ceil, new_y_centre_ceil);
  //(0, x0)
  output.bottomRightCorner(box_dim_y, box_dim_x) = input.bottomRightCorner(box_dim_y, box_dim_x);
 
  box_dim_x = std::min(old_x_centre_floor, new_x_centre_floor);
  box_dim_y = std::min(old_y_centre_floor, new_y_centre_floor);
  //(y0, x0)
  output.topRightCorner(box_dim_y, box_dim_x) = input.topRightCorner(box_dim_y, box_dim_x);
 
  return output;
}

Referenced by re_sample_image().

re_sample_image()

Matrix< t_complex > purify::utilities::re_sample_image	(	const Matrix< t_complex > &	input,
		const t_real &	re_sample_ratio
	)

resamples image size

Definition at line 365 of file utilities.cc.

                                                                                                 {
  const auto ft_plan = operators::fftw_plan::estimate;
  auto fft =
      purify::operators::init_FFT_2d<Vector<t_complex>>(input.rows(), input.cols(), 1., ft_plan);
  Vector<t_complex> ft_grid = Vector<t_complex>::Zero(input.size());
  std::get<0>(fft)(ft_grid, Vector<t_complex>::Map(input.data(), input.size()));
  Matrix<t_complex> const new_ft_grid = utilities::re_sample_ft_grid(
      Matrix<t_complex>::Map(ft_grid.data(), input.rows(), input.cols()), re_sample_ratio);
  Vector<t_complex> output = Vector<t_complex>::Zero(input.size());
  fft = purify::operators::init_FFT_2d<Vector<t_complex>>(new_ft_grid.rows(), new_ft_grid.cols(),
                                                          1., ft_plan);
  std::get<1>(fft)(output, Vector<t_complex>::Map(new_ft_grid.data(), new_ft_grid.size()));
  return Matrix<t_complex>::Map(output.data(), new_ft_grid.rows(), new_ft_grid.cols()) *
         re_sample_ratio;
}

References purify::operators::estimate, and re_sample_ft_grid().

read_ant_positions()

Matrix< t_real > purify::utilities::read_ant_positions

(

const std::string &

pos_name

)

Read in a text file of antenna positions into a matrix [x, y ,z].

Definition at line 77 of file uvw_utilities.cc.

                                                             {
  /*
    Reads an csv file with x, y, z (meters) and returns the vectors as a matrix (x, y, z) cols.
 
    vis_name:: name of input text file containing [x, y, z] (separated by ' ').
  */
  std::ifstream pos_file(pos_name);
  pos_file.precision(13);
  t_int row = 0;
  std::string line;
  // counts size of pos file
  while (std::getline(pos_file, line)) ++row;
  if (row < 1) throw std::runtime_error("No positions in the file: " + pos_name);
  Matrix<t_real> B = Matrix<t_real>::Zero(row, 3);
 
  pos_file.clear();
  pos_file.seekg(0);
  // reads in vis file
  for (row = 0; row < B.rows(); ++row) {
    B(row, 0) = streamtoreal(pos_file);
    B(row, 1) = streamtoreal(pos_file);
    B(row, 2) = streamtoreal(pos_file);
  }
 
  return B;
}

References streamtoreal().

Referenced by read_ant_positions_to_coverage(), and TEST_CASE().

read_ant_positions_to_coverage()

template<class T , class K >

utilities::vis_params purify::utilities::read_ant_positions_to_coverage	(	const std::string &	pos_name,
		const T &	frequencies,
		const K &	times,
		const t_real	theta_ra,
		const t_real	phi_dec,
		const t_real	latitude
	)

Read in a text file of antenna positions into a matrix [x, y ,z] and generate coverage for frequencies

Definition at line 183 of file uvw_utilities.h.

                                                                            {
  return antenna_to_coverage(read_ant_positions(pos_name), frequencies, times, theta_ra, phi_dec,
                             latitude);
}

References antenna_to_coverage(), and read_ant_positions().

Referenced by main(), and TEST_CASE().

read_visibility() [1/3]

utilities::vis_params purify::utilities::read_visibility	(	const std::string &	vis_name,
		const bool	w_term
	)

Reads in visibility file.

Definition at line 236 of file uvw_utilities.cc.

                                                                                  {
#ifdef PURIFY_H5
  if (has_suffix(vis_name, ".h5")) return H5::read_visibility(vis_name, w_term);
#endif
  return read_visibility_csv(vis_name, w_term);
}

References has_suffix(), purify::H5::read_visibility(), and read_visibility_csv().

read_visibility() [2/3]

utilities::vis_params purify::utilities::read_visibility	(	const std::string &	vis_name2,
		const utilities::vis_params &	uv1
	)

Reads in two visibility files.

Definition at line 150 of file uvw_utilities.cc.

                                                                      {
  const bool w_term = not uv1.w.isZero(0);
 
  const auto uv2 = read_visibility(vis_name2, w_term);
  utilities::vis_params uv;
  uv.u = Vector<t_real>::Zero(uv1.size() + uv2.size());
  uv.v = Vector<t_real>::Zero(uv1.size() + uv2.size());
  uv.w = Vector<t_real>::Zero(uv1.size() + uv2.size());
  uv.vis = Vector<t_complex>::Zero(uv1.size() + uv2.size());
  uv.weights = Vector<t_complex>::Zero(uv1.size() + uv2.size());
  uv.u.segment(0, uv1.size()) = uv1.u;
  uv.v.segment(0, uv1.size()) = uv1.v;
  uv.w.segment(0, uv1.size()) = uv1.w;
  uv.vis.segment(0, uv1.size()) = uv1.vis;
  uv.weights.segment(0, uv1.size()) = uv1.weights;
  uv.u.segment(uv1.size(), uv2.size()) = uv2.u;
  uv.v.segment(uv1.size(), uv2.size()) = uv2.v;
  uv.w.segment(uv1.size(), uv2.size()) = uv2.w;
  uv.vis.segment(uv1.size(), uv2.size()) = uv2.vis;
  uv.weights.segment(uv1.size(), uv2.size()) = uv2.weights;
  return uv;
}

References read_visibility(), purify::utilities::vis_params::size(), purify::utilities::vis_params::u, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

read_visibility() [3/3]

utilities::vis_params purify::utilities::read_visibility	(	const std::vector< std::string > &	names,
		const bool	w_term
	)

Read visibility files from name of vector.

Definition at line 144 of file uvw_utilities.cc.

                                                                                          {
  utilities::vis_params output = read_visibility(names.at(0), w_term);
  if (names.size() == 1) return output;
  for (int i = 1; i < names.size(); i++) output = read_visibility(names.at(i), output);
  return output;
}

Referenced by dirty_visibilities(), main(), b_utilities::random_measurements(), purify::read_measurements::read_measurements(), read_visibility(), and TEST_CASE().

read_visibility_csv()

utilities::vis_params purify::utilities::read_visibility_csv	(	const std::string &	vis_name,
		const bool	w_term
	)

Reads in visibility csv file.

Definition at line 181 of file uvw_utilities.cc.

                                                                                      {
  /*
    Reads an csv file with u, v, visibilities and returns the vectors.
 
    vis_name:: name of input text file containing [u, v, real(V), imag(V)] (separated by ' ').
  */
  std::ifstream vis_file(vis_name);
  if (vis_file) {
    PURIFY_LOW_LOG("File {} successfully opened", vis_name);
  } else {
    throw std::runtime_error("Could not open file " + vis_name);
  }
  vis_file.precision(13);
  t_int row = 0;
  std::string line;
  // counts size of vis file
  while (std::getline(vis_file, line)) ++row;
  Vector<t_real> utemp = Vector<t_real>::Zero(row);
  Vector<t_real> vtemp = Vector<t_real>::Zero(row);
  Vector<t_real> wtemp = Vector<t_real>::Zero(row);
  Vector<t_complex> vistemp = Vector<t_complex>::Zero(row);
  Vector<t_complex> weightstemp = Vector<t_complex>::Zero(row);
 
  vis_file.clear();
  vis_file.seekg(0);
  // reads in vis file
  t_real real;
  t_real imag;
  t_real entry;
  for (row = 0; row < vistemp.size(); ++row) {
    utemp(row) = streamtoreal(vis_file);
    vtemp(row) = streamtoreal(vis_file);
    if (w_term) {
      wtemp(row) = streamtoreal(vis_file);
    }
    real = streamtoreal(vis_file);
    imag = streamtoreal(vis_file);
    entry = streamtoreal(vis_file);
    vistemp(row) = t_complex(real, imag);
    weightstemp(row) = 1 / entry;
  }
  utilities::vis_params uv_vis;
  uv_vis.u = utemp;
  uv_vis.v = -vtemp;  // found that a reflection is needed for the orientation of the gridded image
                      // to be correct
  uv_vis.w = wtemp;
  uv_vis.vis = vistemp;
  uv_vis.weights = weightstemp;
  uv_vis.ra = 0;
  uv_vis.dec = 0;
  uv_vis.average_frequency = 0;
 
  return uv_vis;
}

References purify::utilities::vis_params::average_frequency, purify::utilities::vis_params::dec, PURIFY_LOW_LOG, purify::utilities::vis_params::ra, streamtoreal(), purify::utilities::vis_params::u, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

Referenced by read_visibility().

regroup() [1/2]

void purify::utilities::regroup	(	vis_params &	uv_params,
		std::vector< t_int > &	image_index,
		std::vector< t_int > const &	groups_,
		const t_int	max_groups
	)

Definition at line 14 of file mpi_utilities.cc.

                                                                      {
  std::vector<t_int> groups = groups_;
  // Figure out size of each group
  std::map<t_int, t_int> sizes;
  for (auto g = 0; g < max_groups; g++) sizes[g] = 0;
  for (auto const &group : groups) ++sizes[group];
 
  std::map<t_int, t_int> indices, ends;
  auto i = 0;
  for (auto const &item : sizes) {
    indices[item.first] = i;
    i += item.second;
    ends[item.first] = i;
  }
  const auto minmax = std::minmax_element(ends.begin(), ends.end());
  if (std::get<1>(*std::get<0>(minmax)) < 0)
    throw std::runtime_error("segment end " + std::to_string(std::get<1>(*std::get<0>(minmax))) +
                             " less than 0. Not a valid end.");
  if (std::get<1>(*std::get<1>(minmax)) > uv_params.size())
    throw std::runtime_error("segment end " + std::to_string(std::get<1>(*std::get<1>(minmax))) +
                             " larger than data vector " +
                             std::to_string(static_cast<t_int>(uv_params.size())) +
                             ". End not valid.");
 
  auto const expected = [&ends](t_int i) {
    t_int j = 0;
    for (auto const end : ends) {
      if (i < end.second) return j;
      ++j;
    }
    return j;
  };
 
  i = 0;
  while (i < uv_params.u.size()) {
    auto const expected_proc = expected(i);
    if (groups[i] == expected_proc) {
      ++i;
      continue;
    }
    auto &swapper = indices[groups[i]];
    if (groups[swapper] == expected(swapper)) {
      ++swapper;
      continue;
    }
    if (swapper >= uv_params.u.size())
      throw std::runtime_error("regroup (groups, " + std::to_string(groups[swapper]) + ", " +
                               std::to_string(groups[i]) + ") index out of bounds " +
                               std::to_string(i) + " " + std::to_string(swapper) +
                               " >= " + std::to_string(uv_params.u.size()));
    std::swap(groups[i], groups[swapper]);
    std::swap(uv_params.u(i), uv_params.u(swapper));
    std::swap(uv_params.v(i), uv_params.v(swapper));
    std::swap(uv_params.w(i), uv_params.w(swapper));
    std::swap(uv_params.vis(i), uv_params.vis(swapper));
    std::swap(uv_params.weights(i), uv_params.weights(swapper));
    std::swap(image_index[i], image_index[swapper]);
 
    ++swapper;
  }
}

References purify::utilities::vis_params::size(), purify::utilities::vis_params::u, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

regroup() [2/2]

void purify::utilities::regroup	(	vis_params &	uv_params,
		std::vector< t_int > const &	groups_,
		const t_int	max_groups
	)

Definition at line 9 of file mpi_utilities.cc.

                                                                                             {
  std::vector<t_int> image_index(uv_params.size(), 0);
  regroup(uv_params, image_index, groups_, max_groups);
}

References purify::utilities::vis_params::size().

Referenced by regroup_and_all_to_all(), regroup_and_scatter(), and TEST_CASE().

regroup_and_all_to_all() [1/2]

std::tuple<vis_params, std::vector<t_int> > purify::utilities::regroup_and_all_to_all	(	vis_params const &	params,
		const std::vector< t_int > &	image_index,
		std::vector< t_int > const &	groups,
		sopt::mpi::Communicator const &	comm
	)

Definition at line 94 of file mpi_utilities.cc.

                                                                     {
  if (comm.size() == 1) return std::make_tuple(params, image_index);
  vis_params copy = params;
  std::vector<t_int> index_copy(image_index.size());
  std::copy(image_index.begin(), image_index.end(), index_copy.begin());
  regroup(copy, index_copy, groups, comm.size());
 
  std::vector<t_int> sizes(comm.size());
  std::fill(sizes.begin(), sizes.end(), 0);
  for (const t_int &group : groups) {
    if (group >= static_cast<t_int>(comm.size())) {
      throw std::out_of_range("groups should go from 0 to comm.size()-1");
    }
    ++sizes[group];
  }
 
  return std::make_tuple(all_to_all_visibilities(copy, sizes, comm),
                         comm.all_to_allv(index_copy, sizes));
}

References all_to_all_visibilities(), regroup(), and purify::utilities::vis_params::size().

Referenced by regroup_and_all_to_all(), TEST_CASE(), w_stacking(), and w_stacking_with_all_to_all().

regroup_and_all_to_all() [2/2]

vis_params purify::utilities::regroup_and_all_to_all	(	vis_params const &	params,
		std::vector< t_int > const &	groups,
		sopt::mpi::Communicator const &	comm
	)

Definition at line 116 of file mpi_utilities.cc.

                                                                   {
  return std::get<0>(
      regroup_and_all_to_all(params, std::vector<t_int>(params.size(), 0), groups, comm));
}

References regroup_and_all_to_all(), and purify::utilities::vis_params::size().

regroup_and_scatter()

vis_params purify::utilities::regroup_and_scatter	(	vis_params const &	params,
		std::vector< t_int > const &	groups,
		sopt::mpi::Communicator const &	comm
	)

Definition at line 77 of file mpi_utilities.cc.

                                                                {
  if (comm.size() == 1) return params;
  if (comm.rank() != comm.root_id()) return scatter_visibilities(comm);
 
  std::vector<t_int> sizes(comm.size());
  std::fill(sizes.begin(), sizes.end(), 0);
  for (auto const &group : groups) {
    if (group > static_cast<t_int>(comm.size()))
      throw std::out_of_range("groups should go from 0 to comm.size()");
    ++sizes[group];
  }
 
  vis_params copy = params;
  regroup(copy, groups, comm.size());
  return scatter_visibilities(copy, sizes, comm);
}

References regroup(), scatter_visibilities(), and purify::utilities::vis_params::size().

Referenced by dirty_visibilities(), distribute_params(), and TEST_CASE().

scatter_visibilities() [1/2]

vis_params purify::utilities::scatter_visibilities

(

sopt::mpi::Communicator const &

comm

)

Definition at line 156 of file mpi_utilities.cc.

                                                                 {
  if (comm.is_root())
    throw std::runtime_error("The root node should call the *other* scatter_visibilities function");
 
  auto const local_size = comm.scatter_one<t_int>();
  vis_params result;
  result.u = comm.scatterv<decltype(result.u)::Scalar>(local_size);
  result.v = comm.scatterv<decltype(result.v)::Scalar>(local_size);
  result.w = comm.scatterv<decltype(result.w)::Scalar>(local_size);
  result.vis = comm.scatterv<decltype(result.vis)::Scalar>(local_size);
  result.weights = comm.scatterv<decltype(result.weights)::Scalar>(local_size);
  result.units = static_cast<utilities::vis_units>(
      comm.broadcast<std::remove_const<decltype(static_cast<int>(result.units))>::type>());
  result.ra = comm.broadcast<std::remove_const<decltype(result.ra)>::type>();
  result.dec = comm.broadcast<std::remove_const<decltype(result.dec)>::type>();
  result.average_frequency =
      comm.broadcast<std::remove_const<decltype(result.average_frequency)>::type>();
  return result;
}

References purify::utilities::vis_params::u.

scatter_visibilities() [2/2]

vis_params purify::utilities::scatter_visibilities	(	vis_params const &	params,
		std::vector< t_int > const &	sizes,
		sopt::mpi::Communicator const &	comm
	)

Definition at line 137 of file mpi_utilities.cc.

                                                                 {
  if (comm.size() == 1) return params;
  if (not comm.is_root()) return scatter_visibilities(comm);
 
  comm.scatter_one(sizes);
  vis_params result;
  result.u = comm.scatterv(params.u, sizes);
  result.v = comm.scatterv(params.v, sizes);
  result.w = comm.scatterv(params.w, sizes);
  result.vis = comm.scatterv(params.vis, sizes);
  result.weights = comm.scatterv(params.weights, sizes);
  result.units = static_cast<utilities::vis_units>(comm.broadcast(static_cast<int>(params.units)));
  result.ra = comm.broadcast(params.ra);
  result.dec = comm.broadcast(params.dec);
  result.average_frequency = comm.broadcast(params.average_frequency);
  return result;
}

References purify::utilities::vis_params::average_frequency, purify::utilities::vis_params::dec, purify::utilities::vis_params::ra, purify::utilities::vis_params::u, purify::utilities::vis_params::units, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

Referenced by dirty_visibilities(), distribute_params(), regroup_and_scatter(), and TEST_CASE().

set_cell_size() [1/3]

utilities::vis_params purify::utilities::set_cell_size	(	const sopt::mpi::Communicator &	comm,
		utilities::vis_params const &	uv_vis,
		const t_real &	cell_x,
		const t_real &	cell_y
	)

Definition at line 186 of file mpi_utilities.cc.

                                                          {
  if (comm.size() == 1) return utilities::set_cell_size(uv_vis, cell_x, cell_y);
 
  const t_real max_u = comm.all_reduce<t_real>(uv_vis.u.array().cwiseAbs().maxCoeff(), MPI_MAX);
  const t_real max_v = comm.all_reduce<t_real>(uv_vis.v.array().cwiseAbs().maxCoeff(), MPI_MAX);
  return utilities::set_cell_size(uv_vis, max_u, max_v, cell_x, cell_y);
}

References purify::utilities::vis_params::u, and purify::utilities::vis_params::v.

Referenced by set_cell_size(), and TEST_CASE().

set_cell_size() [2/3]

utilities::vis_params purify::utilities::set_cell_size	(	const utilities::vis_params &	uv_vis,
		const t_real &	cell_size_u,
		const t_real &	cell_size_v
	)

Scales visibilities to a given pixel size in arcseconds.

Definition at line 304 of file uvw_utilities.cc.

                                                               {
  const t_real max_u = std::sqrt((uv_vis.u.array() * uv_vis.u.array()).maxCoeff());
  const t_real max_v = std::sqrt((uv_vis.v.array() * uv_vis.v.array()).maxCoeff());
  return set_cell_size(uv_vis, max_u, max_v, cell_size_u, cell_size_v);
}

References set_cell_size(), purify::utilities::vis_params::u, and purify::utilities::vis_params::v.

set_cell_size() [3/3]

utilities::vis_params purify::utilities::set_cell_size	(	const utilities::vis_params &	uv_vis,
		const t_real &	max_u,
		const t_real &	max_v,
		const t_real &	input_cell_size_u,
		const t_real &	input_cell_size_v
	)

Definition at line 261 of file uvw_utilities.cc.

                                                                     {
  /*
    Converts the units of visibilities to units of 2 * pi, while scaling for the size of a pixel
    (cell_size)
 
    uv_vis:: visibilities
    cell_size:: size of a pixel in arcseconds
  */
 
  utilities::vis_params scaled_vis = uv_vis;
  t_real cell_size_u = input_cell_size_u;
  t_real cell_size_v = input_cell_size_v;
  if (cell_size_u == 0 and cell_size_v == 0) {
    cell_size_u = (180 * 3600) / max_u / constant::pi / 3;  // Calculate cell size if not given one
 
    cell_size_v = (180 * 3600) / max_v / constant::pi / 3;  // Calculate cell size if not given one
    // PURIFY_MEDIUM_LOG("PSF has a FWHM of {} by {} arcseconds", cell_size_u * 3, cell_size_v * 3);
  }
  if (cell_size_v == 0) {
    cell_size_v = cell_size_u;
  }
 
  PURIFY_MEDIUM_LOG("Using a pixel size of {} by {} arcseconds", cell_size_u, cell_size_v);
  t_real scale_factor_u = 1;
  t_real scale_factor_v = 1;
  if (uv_vis.units == utilities::vis_units::lambda) {
    scale_factor_u = 180 * 3600 / cell_size_u / constant::pi;
    scale_factor_v = 180 * 3600 / cell_size_v / constant::pi;
    scaled_vis.w = uv_vis.w;
  }
  if (uv_vis.units == utilities::vis_units::radians) {
    scale_factor_u = 180 * 3600 / constant::pi;
    scale_factor_v = 180 * 3600 / constant::pi;
    scaled_vis.w = uv_vis.w;
  }
  scaled_vis.u = uv_vis.u / scale_factor_u * 2 * constant::pi;
  scaled_vis.v = uv_vis.v / scale_factor_v * 2 * constant::pi;
 
  scaled_vis.units = utilities::vis_units::radians;
  return scaled_vis;
}

References lambda, purify::constant::pi, PURIFY_MEDIUM_LOG, radians, purify::utilities::vis_params::u, purify::utilities::vis_params::units, purify::utilities::vis_params::v, and purify::utilities::vis_params::w.

SNR_to_standard_deviation()

t_real purify::utilities::SNR_to_standard_deviation	(	const Vector< t_complex > &	y0,
		const t_real &	SNR
	)

Converts SNR to RMS noise.

Definition at line 101 of file utilities.cc.

                                                                                 {
  /*
  Returns value of noise rms given a measurement vector and signal to noise ratio
  y0:: complex valued vector before noise added
  SNR:: signal to noise ratio
 
  This calculation follows Carrillo et al. (2014), PURIFY a new approach to radio interferometric
  imaging but we have assumed that rms noise is for the real and imaginary parts
  */
  return y0.stableNorm() / std::sqrt(2 * y0.size()) * std::pow(10.0, -(SNR / 20.0));
}

Referenced by b_utilities::dirty_measurements(), dirty_visibilities(), getInputData(), and main().

sort_by_w()

vis_params purify::utilities::sort_by_w

(

const vis_params &

uv_data

)

Sort visibilities to be from w_max to w_min.

Definition at line 5 of file wproj_utilities.cc.

                                                {
  vis_params output = uv_data;
  std::vector<t_uint> ind;
  for (t_uint i = 0; i < uv_data.size(); i++) ind.push_back(i);
  std::sort(ind.begin(), ind.end(),
            [&](const t_uint a, const t_uint b) { return uv_data.w(a) < uv_data.w(b); });
  for (t_uint i = 0; i < uv_data.size(); i++) {
    output.u(i) = uv_data.u(ind.at(i));
    output.v(i) = uv_data.v(ind.at(i));
    output.w(i) = uv_data.w(ind.at(i));
    output.weights(i) = uv_data.weights(ind.at(i));
    output.vis(i) = uv_data.vis(ind.at(i));
  }
  return output;
}

References purify::utilities::vis_params::size(), purify::utilities::vis_params::u, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

sparse_multiply_matrix()

template<class T0 , class T1 >

std::enable_if<std::is_same<typename T0::Scalar, typename T1::Scalar>::value and T0::IsRowMajor, Vector<typename T0::Scalar> >::type purify::utilities::sparse_multiply_matrix	(	const Eigen::SparseMatrixBase< T0 > &	M,
		const Eigen::MatrixBase< T1 > &	x
	)

Parallel multiplication with a sparse matrix and vector.

Definition at line 67 of file utilities.h.

                                                                                         {
  assert(M.cols() == x.size());
  Vector<typename T0::Scalar> y = Vector<typename T0::Scalar>::Zero(M.rows());
  auto const &derived = M.derived();
// parallel sparse matrix multiplication with vector.
#pragma omp parallel for
  // #pragma omp simd
  for (t_int k = 0; k < M.outerSize(); ++k)
    for (typename Sparse<typename T0::Scalar>::InnerIterator it(derived, k); it; ++it)
      y(k) += it.value() * x(it.index());
  return y;
}

Referenced by purify::operators::init_gridding_matrix_2d(), and TEST_CASE().

standard_deviation()

template<class K >

t_real purify::utilities::standard_deviation

(

const K

x

)

Calculates the standard deviation of a vector.

Definition at line 34 of file utilities.h.

                                     {
  // calculate standard deviation of vector x
  return std::sqrt(variance(x));
}

References variance().

Referenced by add_noise(), calculate_l2_radius(), generate_antennas(), and random_sample_density().

step_size()

template<class T >

t_real purify::utilities::step_size	(	T const &	vis,
		const std::shared_ptr< sopt::LinearTransform< T > const > &	measurements,
		const std::shared_ptr< sopt::LinearTransform< T > const > &	wavelets,
		const t_uint	sara_size
	)

Calculate step size using MPI (does not include factor of 1e-3)

Parameters

[in]	vis	Vector of measurement data
[in]	measurements	Shared pointer to measurement linear operator
[in]	wavelets	Shared pointer to SARA wavelet linear operator
[in]	sara_size	Size of sara dictionary of MPI node.

Please use this function to calculate the step size for PADMM, Forward Backward, etc. Especially when using MPI.

Note

There is an issue with using the adjoint wavelet transform of the SARA basis. When there are more nodes than wavelets, there will be nodes with an empty Vector<T> of wavelet coefficients. This can cause some functions to break, like .maxCoeff().

Definition at line 79 of file mpi_utilities.h.

                                         {
  // measurement operator may use different number of nodes than wavelet operator
  // so needs to be done separately
  const T dimage = measurements->adjoint() * vis;
  return (sara_size > 0) ? (wavelets->adjoint() * dimage).cwiseAbs().maxCoeff() : 0.;
};

References measurements.

Referenced by purify::factory::fb_factory(), purify::factory::padmm_factory(), padmm_factory(), and purify::factory::primaldual_factory().

streamtoreal()

t_real purify::utilities::streamtoreal

(

std::ifstream &

stream

)

Reading reals from visibility file (including nan's and inf's)

Definition at line 175 of file uvw_utilities.cc.

                                         {
  std::string input;
  stream >> input;
  return std::stod(input);
}

Referenced by read_ant_positions(), and read_visibility_csv().

sub2ind()

t_int purify::utilities::sub2ind	(	const t_int &	row,
		const t_int &	col,
		const t_int &	rows,
		const t_int &	cols
	)

Converts from subscript to index for matrix.

Definition at line 12 of file utilities.cc.

                                                                                        {
  /*
    Converts (row, column) of a matrix to a single index. This does the same as the matlab funciton
    sub2ind, converts subscript to index.
    Though order of cols and rows is probably transposed.
 
    row:: row of matrix (y)
    col:: column of matrix (x)
    cols:: number of columns for matrix
    rows:: number of rows for matrix
   */
  return row * cols + col;
}

Referenced by purify::details::init_gridding_matrix_2d(), purify::operators::init_on_the_fly_gridding_matrix_2d(), purify::operators::init_zero_padding_2d(), purify::wproj_utilities::row_wise_convolution(), and purify::wproj_utilities::row_wise_sparse_convolution().

uv_scale()

utilities::vis_params purify::utilities::uv_scale	(	const utilities::vis_params &	uv_vis,
		const t_int &	sizex,
		const t_int &	sizey
	)

scales the visibilities to units of pixels

Definition at line 311 of file uvw_utilities.cc.

                                                   {
  /*
    scales the uv coordinates from being in units of 2 * pi to units of pixels.
  */
  utilities::vis_params scaled_vis;
  scaled_vis.u = uv_vis.u / (2 * constant::pi) * sizex;
  scaled_vis.v = uv_vis.v / (2 * constant::pi) * sizey;
  scaled_vis.vis = uv_vis.vis;
  scaled_vis.weights = uv_vis.weights;
  for (t_int i = 0; i < uv_vis.u.size(); ++i) {
    // scaled_vis.u(i) = utilities::mod(scaled_vis.u(i), sizex);
    // scaled_vis.v(i) = utilities::mod(scaled_vis.v(i), sizey);
  }
  scaled_vis.w = uv_vis.w;
  scaled_vis.units = utilities::vis_units::pixels;
  scaled_vis.ra = uv_vis.ra;
  scaled_vis.dec = uv_vis.dec;
  scaled_vis.average_frequency = uv_vis.average_frequency;
  return scaled_vis;
}

References purify::utilities::vis_params::average_frequency, purify::utilities::vis_params::dec, purify::constant::pi, pixels, purify::utilities::vis_params::ra, purify::utilities::vis_params::u, purify::utilities::vis_params::units, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

Referenced by TEST_CASE().

variance()

template<class K >

t_real purify::utilities::variance

(

const K

x

)

Calculate variance of vector.

Definition at line 27 of file utilities.h.

                           {
  // calculate variance of vector x
  const Matrix<typename K::Scalar> q = (x.array() - x.array().mean()).matrix();
  return std::real((q.adjoint() * q)(0) / static_cast<t_real>(q.size() - 1));
}

Referenced by standard_deviation(), and TEST_CASE().

w_lambert() [1/2]

template<typename T >

t_real purify::utilities::w_lambert	(	const T &	x,
		const t_real &	relative_difference = 1e-8
	)

Definition at line 292 of file wproj_utilities.h.

                                                                       {
  if (std::is_integral<T>::value) return w_lambert<t_real>(x, relative_difference);
  T z = T(x);
  T diff = T(1e4 * relative_difference);
  if (x < -1 / std::exp(1)) throw std::runtime_error("Out of bounds for W lambert function!");
  while (std::abs(diff) > relative_difference) {
    const T f = z * std::exp(z) - x;
    const T w = z - f / (std::exp(z) * (z + 1) - f * (z + 2) / (2 * z + 2));
    diff = w - z;
    z = w;
  }
  return z;
}

w_lambert() [2/2]

template<typename T >

t_real purify::utilities::w_lambert	(	T	x,
		const t_real &	relative_difference
	)

Calculate W Lambert function.

w_stacking()

utilities::vis_params purify::utilities::w_stacking	(	utilities::vis_params const &	params,
		sopt::mpi::Communicator const &	comm,
		const t_int	iters,
		const std::function< t_real(t_real)> &	cost,
		const t_real	k_means_rel_diff
	)

Definition at line 195 of file mpi_utilities.cc.

                                                                {
  const std::vector<t_int> image_index = std::get<0>(
      distribute::kmeans_algo(params.w, comm.size(), iters, comm, cost, k_means_rel_diff));
  return utilities::regroup_and_all_to_all(params, image_index, comm);
}

References purify::distribute::kmeans_algo(), regroup_and_all_to_all(), and purify::utilities::vis_params::w.

Referenced by purify::operators::base_degrid_operator_2d(), getInputData(), purify::measurementoperator::init_degrid_operator_2d(), main(), and DegridOperatorCtorFixturePar::SetUp().

w_stacking_with_all_to_all()

std::tuple<utilities::vis_params, std::vector<t_int>, std::vector<t_real> > purify::utilities::w_stacking_with_all_to_all	(	utilities::vis_params const &	params,
		const t_real	du,
		const t_int	min_support,
		const t_int	max_support,
		sopt::mpi::Communicator const &	comm,
		const t_int	iters,
		const t_real	fill_relaxation,
		const std::function< t_real(t_real)> &	cost,
		const t_real	k_means_rel_diff
	)

Definition at line 204 of file mpi_utilities.cc.

                                                          {
  const auto kmeans =
      distribute::kmeans_algo(params.w, comm.size(), iters, comm, cost, k_means_rel_diff);
  const std::vector<t_real> &w_stacks = std::get<1>(kmeans);
  const std::vector<t_int> groups = distribute::w_support(
      params.w, std::get<0>(kmeans), w_stacks, du, min_support, max_support, fill_relaxation, comm);
  utilities::vis_params outdata;
  std::vector<t_int> image_index;
  std::tie(outdata, image_index) =
      utilities::regroup_and_all_to_all(params, std::get<0>(kmeans), groups, comm);
  return std::tuple<utilities::vis_params, std::vector<t_int>, std::vector<t_real>>(
      outdata, image_index, w_stacks);
}

References purify::distribute::kmeans_algo(), regroup_and_all_to_all(), purify::utilities::vis_params::w, and purify::widefield::w_support().

Referenced by getInputData(), main(), DegridOperatorCtorFixturePar::SetUp(), and TEST_CASE().

write_visibility()

void purify::utilities::write_visibility	(	const utilities::vis_params &	uv_vis,
		const std::string &	file_name,
		const bool	w_term
	)

Writes visibilities to txt.

Definition at line 243 of file uvw_utilities.cc.

                                         {
  /*
          writes visibilities to output text file (currently ignores w-component)
          uv_vis:: input uv data
          file_name:: name of output text file
  */
  std::ofstream out(file_name);
  out.precision(13);
  for (t_int i = 0; i < uv_vis.u.size(); ++i) {
    out << uv_vis.u(i) << " " << -uv_vis.v(i) << " ";
    if (w_term) out << uv_vis.w(i) << " ";
    out << std::real(uv_vis.vis(i)) << " " << std::imag(uv_vis.vis(i)) << " "
        << 1. / std::real(uv_vis.weights(i)) << std::endl;
  }
  out.close();
}

References purify::utilities::vis_params::u, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

Referenced by main(), padmm(), b_utilities::random_measurements(), and TEST_CASE().