purify::wproj_utilities Namespace Reference

Namespaces
	expansions

Functions
std::vector< t_uint >	w_rows (const Sparse< t_complex > &w_degrider)
std::tuple< t_real, t_real >	fov_cosines (t_real const &cell_x, t_real const &cell_y, t_uint const &x_size, t_uint const &y_size)
	Work out max L and M directional cosines from image parameters. More...
Matrix< t_complex >	generate_dde (const std::function< t_complex(t_real, t_real)> &dde, const t_real &cell_x, const t_real &cell_y, const t_uint &x_size, const t_uint &y_size)
	Generate image of DDE for A or W projection. More...
Matrix< t_complex >	generate_chirp (const std::function< t_complex(t_real, t_real)> &dde, const t_real &w_rate, const t_real &cell_x, const t_real &cell_y, const t_uint &x_size, const t_uint &y_size)
Matrix< t_complex >	generate_chirp (const t_real &w_rate, const t_real &cell_x, const t_real &cell_y, const t_uint &x_size, const t_uint &y_size)
	Generates image of chirp. More...
Sparse< t_complex >	create_chirp_row (const t_real &w_rate, const t_real &cell_x, const t_real &cell_y, const t_uint &ftsizev, const t_uint &ftsizeu, const t_real &energy_fraction, const sopt::OperatorFunction< Vector< t_complex >> &fftop)
	Generates row of chirp matrix from image of chirp. More...
Sparse< t_complex >	create_chirp_row (const Vector< t_complex > &chirp_image, const t_real &energy_fraction, const sopt::OperatorFunction< Vector< t_complex >> &fftop)
Sparse< t_complex >	wprojection_matrix (const Sparse< t_complex > &G, const t_uint &x_size, const t_uint &y_size, const Vector< t_real > &w_components, const t_real &cell_x, const t_real &cell_y, const t_real &energy_fraction_chirp, const t_real &energy_fraction_wproj, const expansions::series series=expansions::series::none, const t_uint order=1, const t_real &interpolation_error=1e-2)
	Produce Gridding matrix convovled with chirp matrix for wprojection. More...
t_real	snr_metric (const Image< t_real > &model, const Image< t_real > &solution)
	SNR calculation. More...
t_real	mr_metric (const Image< t_real > &model, const Image< t_real > &solution)
	MR calculation. More...
Sparse< t_complex >	generate_vect (const t_uint &x_size, const t_uint &y_size, const t_real &sigma, const t_real &mean)
	Genereates an image of a Gaussian as a sparse matrice. More...
template<typename T >
Sparse< t_complex >	convert_sparse (const Eigen::MatrixBase< T > &M, const t_real &threshold=0.)
	Convert from dense to sparse. More...
template<typename T >
t_real	sparsify_row_thres (const Eigen::SparseMatrixBase< T > &row, const t_real &energy)
	Returns threshold to keep a fraction of energy in the sparse row. More...
template<typename T >
t_real	sparsify_row_dense_thres (const Eigen::MatrixBase< T > &row, const t_real &energy)
	Returns threshold to keep a fraction of energy in the dense row. More...
template<class T >
Sparse< t_complex >	row_wise_convolution (const Sparse< t_complex > &Grid_, const Sparse< T > &chirp_, const t_uint &x_size, const t_uint &y_size)
	Perform convolution with gridding matrix row and chirp matrix row. More...
t_real	sparsity_sp (const Sparse< t_complex > &Gmat)
	return fraction of non zero values from sparse matrix More...
t_real	sparsity_im (const Image< t_complex > &Cmat)
	return faction of non zero values from matrix More...
t_real	upsample_ratio_sim (const utilities::vis_params &uv_vis, const t_real &L, const t_real &M, const t_int &x_size, const t_int &y_size, const t_int &multipleOf)
	Calculate upsample ratio from bandwidth (only needed for simulations) More...
template<class T >
Sparse< t_complex >	row_wise_sparse_convolution (const Sparse< t_complex > &Grid_, const Sparse< T > &chirp_, const t_uint &x_size, const t_uint &y_size)

Function Documentation

convert_sparse()

template<typename T >

Sparse< t_complex > purify::wproj_utilities::convert_sparse	(	const Eigen::MatrixBase< T > &	M,
		const t_real &	threshold = 0.
	)

Convert from dense to sparse.

Definition at line 90 of file wproj_utilities.h.

                                                                                       {
  const Sparse<t_complex> M_sparse =
      M.sparseView(std::max(M.cwiseAbs().maxCoeff() * 1e-10, threshold), 1);
  Sparse<t_complex> out(M.rows(), M.cols());
  out.reserve(M_sparse.nonZeros());
  for (t_int k = 0; k < M_sparse.outerSize(); ++k)
    for (Sparse<t_complex>::InnerIterator it(M_sparse, k); it; ++it) {
      out.coeffRef(it.row(), it.col()) = it.value();
    }
  return out;
}

Referenced by create_chirp_row(), and row_wise_convolution().

create_chirp_row() [1/2]

Sparse< t_complex > purify::wproj_utilities::create_chirp_row	(	const t_real &	w_rate,
		const t_real &	cell_x,
		const t_real &	cell_y,
		const t_uint &	ftsizev,
		const t_uint &	ftsizeu,
		const t_real &	energy_fraction,
		const sopt::OperatorFunction< Vector< t_complex >> &	fftop
	)

Generates row of chirp matrix from image of chirp.

Definition at line 110 of file wproj_utilities.cc.

                                                                                         {
  const Matrix<t_complex> chirp =
      wproj_utilities::generate_chirp(w_rate, cell_x, cell_y, ftsizeu, ftsizev);
  return create_chirp_row(Vector<t_complex>::Map(chirp.data(), chirp.size()), energy_fraction,
                          fftop);
}

References generate_chirp().

Referenced by wprojection_matrix().

create_chirp_row() [2/2]

Sparse< t_complex > purify::wproj_utilities::create_chirp_row	(	const Vector< t_complex > &	chirp_image,
		const t_real &	energy_fraction,
		const sopt::OperatorFunction< Vector< t_complex >> &	fftop
	)

Definition at line 119 of file wproj_utilities.cc.

                                                                                         {
  Vector<t_complex> chirp;
#pragma omp critical(fft)
  fftop(chirp, chirp_image);
  chirp *= std::sqrt(chirp.size());
  const t_real thres = wproj_utilities::sparsify_row_dense_thres(chirp, energy_fraction);
  assert(chirp.norm() > 0);
  const t_real row_max = chirp.cwiseAbs().maxCoeff();
  return convert_sparse(chirp, thres).transpose();
}

References convert_sparse(), and sparsify_row_dense_thres().

fov_cosines()

std::tuple< t_real, t_real > purify::wproj_utilities::fov_cosines	(	t_real const &	cell_x,
		t_real const &	cell_y,
		t_uint const &	x_size,
		t_uint const &	y_size
	)

Work out max L and M directional cosines from image parameters.

Definition at line 37 of file wproj_utilities.cc.

                                                                                   {
  // celly:: size of y pixel in arcseconds
  // cellx:: size of x pixel in arcseconds
  // x_size:: number of pixels along x-axis
  // y_size:: number of pixels along y-axis
  const t_real pi = constant::pi;
  const t_real theta_FoV_L = cell_x * x_size;
  const t_real theta_FoV_M = cell_y * y_size;
  const t_real L = 2 * std::sin(pi / 180. * theta_FoV_L / (60. * 60.) * 0.5);
  const t_real M = 2 * std::sin(pi / 180. * theta_FoV_M / (60. * 60.) * 0.5);
  return std::make_tuple(L, M);
}

References purify::constant::pi.

Referenced by generate_dde().

generate_chirp() [1/2]

Matrix< t_complex > purify::wproj_utilities::generate_chirp	(	const std::function< t_complex(t_real, t_real)> &	dde,
		const t_real &	w_rate,
		const t_real &	cell_x,
		const t_real &	cell_y,
		const t_uint &	x_size,
		const t_uint &	y_size
	)

Definition at line 76 of file wproj_utilities.cc.

                                                                             {
  // return chirp image fourier transform for w component
  // w:: w-term in units of lambda
  // celly:: size of y pixel in arcseconds
  // cellx:: size of x pixel in arcseconds
  // x_size:: number of pixels along x-axis
  // y_size:: number of pixels along y-axis
  const t_real nz = y_size * x_size;
  const t_complex I(0, 1);
  const t_real pi = constant::pi;
  const auto chirp = [=](const t_real &y, const t_real &x) {
    return dde(y, x) * (std::exp(-2 * pi * I * w_rate * (std::sqrt(1 - x * x - y * y) - 1))) / nz;
  };
  auto primary_beam = [=](const t_real &x, const t_real &y) {
    return 1. / std::sqrt(1 - x * x - y * y);
  };
  const auto chirp_approx = [=](const t_real &y, const t_real &x) {
    return dde(y, x) * (std::exp(2 * pi * I * w_rate * (x * x + y * y))) / nz;
  };
  return generate_dde(chirp, cell_x, cell_y, x_size, y_size);
}

References generate_dde(), purify::I, and purify::constant::pi.

Referenced by create_chirp_row(), and generate_chirp().

generate_chirp() [2/2]

Matrix< t_complex > purify::wproj_utilities::generate_chirp	(	const t_real &	w_rate,
		const t_real &	cell_x,
		const t_real &	cell_y,
		const t_uint &	x_size,
		const t_uint &	y_size
	)

Generates image of chirp.

Definition at line 99 of file wproj_utilities.cc.

                                                                             {
  // return chirp image fourier transform for w component
  // w:: w-term in units of lambda
  // celly:: size of y pixel in arcseconds
  // cellx:: size of x pixel in arcseconds
  // x_size:: number of pixels along x-axis
  // y_size:: number of pixels along y-axis
  return generate_chirp([](const t_real &, const t_real &) { return 1.; }, w_rate, cell_x, cell_y,
                        x_size, y_size);
}

References generate_chirp().

generate_dde()

Matrix< t_complex > purify::wproj_utilities::generate_dde	(	const std::function< t_complex(t_real, t_real)> &	dde,
		const t_real &	cell_x,
		const t_real &	cell_y,
		const t_uint &	x_size,
		const t_uint &	y_size
	)

Generate image of DDE for A or W projection.

Definition at line 51 of file wproj_utilities.cc.

                                                     {
  // celly:: size of y pixel in arcseconds
  // cellx:: size of x pixel in arcseconds
  // x_size:: number of pixels along x-axis
  // y_size:: number of pixels along y-axis
  auto const LM = fov_cosines(cell_x, cell_y, x_size, y_size);
 
  const t_real L = std::get<0>(LM);
  const t_real M = std::get<1>(LM);
 
  const t_real delt_x = L / x_size;
  const t_real delt_y = M / y_size;
  Image<t_complex> chirp(y_size, x_size);
 
  for (t_int l = 0; l < x_size; ++l)
    for (t_int m = 0; m < y_size; ++m) {
      const t_real x = (l + 0.5 - x_size * 0.5) * delt_x;
      const t_real y = (m + 0.5 - y_size * 0.5) * delt_y;
      chirp(l, m) = dde(y, x);
    }
 
  return chirp;
}

References fov_cosines().

Referenced by generate_chirp().

generate_vect()

Sparse< t_complex > purify::wproj_utilities::generate_vect	(	const t_uint &	x_size,
		const t_uint &	y_size,
		const t_real &	sigma,
		const t_real &	mean
	)

Genereates an image of a Gaussian as a sparse matrice.

Definition at line 213 of file wproj_utilities.cc.

                                                    {
  // t_real sigma = 3;
  // const t_real pi = constant::pi;
  t_int W = x_size;
  Matrix<t_complex> chirp = Image<t_complex>::Zero(y_size, x_size);
  // t_real mean = 0;
  t_complex I(0, 1);
  t_real step = W / 2;
 
  // t_complex sum = 0.0; // For accumulating the kernel values
  for (int x = 0; x < W; ++x) {
    for (int y = 0; y < W; ++y) {
      t_int xx = x - step;
      xx = xx * xx;
      t_int yy = y - step;
      yy = yy * yy;
      t_complex val =
          exp(-0.5 * (xx / (sigma * sigma) +
                      yy / (sigma * sigma)));  //* std::exp(- 2 * pi * I * (x * 0.5 + y * 0.5));
      if (std::abs(val) > 1e-5)
        chirp(x, y) = val;
      else
        chirp(x, y) = 0.0;
    }
  }
 
  chirp.resize(x_size * y_size, 1);
  Sparse<t_complex> chirp_ = chirp.sparseView(1e-5, 1);
  return chirp_;
}

References purify::I.

mr_metric()

t_real purify::wproj_utilities::mr_metric	(	const Image< t_real > &	model,
		const Image< t_real > &	solution
	)

MR calculation.

Definition at line 199 of file wproj_utilities.cc.

                                                                            {
  /*
     Returns SNR of the estimated model image
     */
  t_int Npix = model.rows() * model.cols();
  Image<t_real> model_ = model.array() + 1e-10;
  Image<t_real> solution_ = solution.array() + 1e-10;
  Image<t_real> model_sol = model_.matrix().cwiseQuotient(solution_.matrix());
  Image<t_real> sol_model = solution_.matrix().cwiseQuotient(model_.matrix());
  Image<t_real> min_ratio = sol_model.matrix().cwiseMin(model_sol.matrix());
  t_real val = min_ratio.array().sum() / Npix;
  return val;
}

row_wise_convolution()

template<class T >

Sparse< t_complex > purify::wproj_utilities::row_wise_convolution	(	const Sparse< t_complex > &	Grid_,
		const Sparse< T > &	chirp_,
		const t_uint &	x_size,
		const t_uint &	y_size
	)

Perform convolution with gridding matrix row and chirp matrix row.

Definition at line 187 of file wproj_utilities.h.

                                                                                   {
  assert((Grid_.nonZeros() > 0) and (chirp_.nonZeros() > 0));
  Matrix<t_complex> output_row = Matrix<t_complex>::Zero(1, x_size * y_size);
  const t_int x_sizec = std::floor(x_size * 0.5);
  const t_int y_sizec = std::floor(y_size * 0.5);
  for (t_int k = 0; k < output_row.size(); k++) {
    // shifting indicies to workout linear convolution
    t_int i_shift_W = 0;
    t_int j_shift_W = 0;
    std::tie(j_shift_W, i_shift_W) = utilities::ind2sub(k, y_size, x_size);
    const t_int i_W = utilities::mod(i_shift_W + x_sizec, x_size);
    const t_int j_W = utilities::mod(j_shift_W + y_sizec, y_size);
    for (Sparse<t_complex>::InnerIterator it(Grid_, 0); it; ++it) {
      t_int i_shift_G = 0;
      t_int j_shift_G = 0;
      std::tie(j_shift_G, i_shift_G) = utilities::ind2sub(it.index(), y_size, x_size);
      const t_int i_G = utilities::mod(i_shift_G + x_sizec, x_size);
      const t_int j_G = utilities::mod(j_shift_G + y_sizec, y_size);
      const t_int i_C = i_W - i_G;
      const t_int j_C = j_W - j_G;
      // Checking that convolution result is going to be within bounds
      if (-x_sizec <= i_C and i_C < std::ceil(x_size * 0.5) and -y_sizec <= j_C and
          j_C < std::ceil(y_size * 0.5)) {
        // FFTshift of result to go into gridding matrix
        const t_int i_shift_C = utilities::mod(i_C, x_size);
        const t_int j_shift_C = utilities::mod(j_C, y_size);
        const t_int pos = utilities::sub2ind(j_shift_C, i_shift_C, y_size, x_size);
        assert(0 <= pos and pos < x_size * y_size);
        output_row(0, k) += it.value() * chirp_.coeff(0, pos);
      }
    }
  }
  return convert_sparse(output_row);
}

References convert_sparse(), purify::utilities::ind2sub(), purify::utilities::mod(), and purify::utilities::sub2ind().

row_wise_sparse_convolution()

template<class T >

Sparse<t_complex> purify::wproj_utilities::row_wise_sparse_convolution	(	const Sparse< t_complex > &	Grid_,
		const Sparse< T > &	chirp_,
		const t_uint &	x_size,
		const t_uint &	y_size
	)

Definition at line 223 of file wproj_utilities.h.

                                                                    {
  assert((Grid_.nonZeros() > 0) and (chirp_.nonZeros() > 0));
  const t_int x_sizec = std::floor(x_size * 0.5);
  const t_int y_sizec = std::floor(y_size * 0.5);
  std::set<t_int> non_zero_W;
  for (Sparse<t_complex>::InnerIterator jt(chirp_, 0); jt; ++jt) {
    // shifting indicies to workout linear convolution
    t_int i_shift_C = 0;
    t_int j_shift_C = 0;
    std::tie(j_shift_C, i_shift_C) = utilities::ind2sub(jt.index(), y_size, x_size);
    const t_int i_C = utilities::mod(i_shift_C + x_sizec, x_size) - x_sizec;
    const t_int j_C = utilities::mod(j_shift_C + y_sizec, y_size) - y_sizec;
    for (Sparse<t_complex>::InnerIterator it(Grid_, 0); it; ++it) {
      t_int i_shift_G = 0;
      t_int j_shift_G = 0;
      std::tie(j_shift_G, i_shift_G) = utilities::ind2sub(it.index(), y_size, x_size);
      const t_int i_G = utilities::mod(i_shift_G + x_sizec, x_size) - x_sizec;
      const t_int j_G = utilities::mod(j_shift_G + y_sizec, y_size) - y_sizec;
      const t_int i_W = i_G + i_C;
      const t_int j_W = j_G + j_C;
      if (-x_sizec <= i_W and i_W < std::ceil(x_size * 0.5) and -y_sizec <= j_W and
          j_W < std::ceil(y_size * 0.5)) {
        const t_int i_shift_W = utilities::mod(i_W, x_size);
        const t_int j_shift_W = utilities::mod(j_W, y_size);
        const t_real pos = utilities::sub2ind(j_shift_W, i_shift_W, y_size, x_size);
        non_zero_W.insert(pos);
      }
    }
  }
  if (non_zero_W.size() < std::min(Grid_.nonZeros(), chirp_.nonZeros()))
    throw std::runtime_error("Gridding kernel or chirp kernel is zero.");
  Sparse<t_complex> output_row(1, x_size * y_size);
  output_row.reserve(non_zero_W.size());
 
  for (t_int k = 0; k < output_row.size(); k++) {
    // shifting indicies to workout linear convolution
    t_int i_shift_W = 0;
    t_int j_shift_W = 0;
    std::tie(j_shift_W, i_shift_W) = utilities::ind2sub(k, y_size, x_size);
    const t_int i_W = utilities::mod(i_shift_W + x_sizec, x_size);
    const t_int j_W = utilities::mod(j_shift_W + y_sizec, y_size);
    for (Sparse<t_complex>::InnerIterator it(Grid_, 0); it; ++it) {
      t_int i_shift_G = 0;
      t_int j_shift_G = 0;
      std::tie(j_shift_G, i_shift_G) = utilities::ind2sub(it.index(), y_size, x_size);
      const t_int i_G = utilities::mod(i_shift_G + x_sizec, x_size);
      const t_int j_G = utilities::mod(j_shift_G + y_sizec, y_size);
      const t_int i_C = i_W - i_G;
      const t_int j_C = j_W - j_G;
      // Checking that convolution result is going to be within bounds
      if (-x_sizec <= i_C and i_C < std::ceil(x_size * 0.5) and -y_sizec <= j_C and
          j_C < std::ceil(y_size * 0.5)) {
        // FFTshift of result to go into gridding matrix
        const t_int i_shift_C = utilities::mod(i_C, x_size);
        const t_int j_shift_C = utilities::mod(j_C, y_size);
        const t_int pos = utilities::sub2ind(j_shift_C, i_shift_C, y_size, x_size);
        assert(0 <= pos and pos < x_size * y_size);
        output_row.coeffRef(0, k) += it.value() * chirp_.coeff(0, pos);
      }
    }
  }
  return output_row;
}

References purify::utilities::ind2sub(), purify::utilities::mod(), and purify::utilities::sub2ind().

Referenced by wprojection_matrix().

snr_metric()

t_real purify::wproj_utilities::snr_metric	(	const Image< t_real > &	model,
		const Image< t_real > &	solution
	)

SNR calculation.

Definition at line 190 of file wproj_utilities.cc.

                                                                             {
  /*
     Returns SNR of the estimated model image
     */
  t_real nm = model.matrix().norm();
  t_real ndiff = (model - solution).matrix().norm();
  t_real val = 20 * std::log10(nm / ndiff);
  return val;
}

sparsify_row_dense_thres()

template<typename T >

t_real purify::wproj_utilities::sparsify_row_dense_thres	(	const Eigen::MatrixBase< T > &	row,
		const t_real &	energy
	)

Returns threshold to keep a fraction of energy in the dense row.

Definition at line 147 of file wproj_utilities.h.

                                                                                     {
  /*
     Takes in a row of G and returns indexes of coeff to keep in the row sparse version
     energy:: how much energy - in l2 sens - to keep after hard-thresholding
     */
 
  if (energy >= 1) return 0;
 
  const Vector<t_real> row_abs = row.cwiseAbs();
  const t_real abs_row_max = row.cwiseAbs().maxCoeff();
  const t_real row_total_energy = row_abs.cwiseProduct(row_abs).sum();
  const t_real target_threshold_energy_sum = row_total_energy * energy;
  const t_int niters = 200;
  t_real lower = 0;
  t_real upper = abs_row_max;
  t_real current_threshold = upper * 0.5;
  t_real old_threshold = upper;
  bool converged = false;
  for (t_int i = 0; i < niters; i++) {
    t_real threshold_energy_sum = 0;
    for (t_int j = 0; j < row_abs.size(); j++) {
      if (row_abs(j) > current_threshold) threshold_energy_sum += row_abs(j) * row_abs(j);
    }
    if ((std::abs(current_threshold - old_threshold) / std::abs(old_threshold) < 1e-3) and
        ((threshold_energy_sum >= target_threshold_energy_sum) and
         (((threshold_energy_sum / target_threshold_energy_sum) - 1) < 0.01))) {
      converged = true;
      break;
    }
    if (threshold_energy_sum < target_threshold_energy_sum)
      upper = current_threshold;
    else
      lower = current_threshold;
    old_threshold = current_threshold;
    current_threshold = (lower + upper) * 0.5;
  }
  if (!converged) current_threshold = lower;
  return current_threshold;
}

Referenced by create_chirp_row().

sparsify_row_thres()

template<typename T >

t_real purify::wproj_utilities::sparsify_row_thres	(	const Eigen::SparseMatrixBase< T > &	row,
		const t_real &	energy
	)

Returns threshold to keep a fraction of energy in the sparse row.

Definition at line 102 of file wproj_utilities.h.

                                                                                     {
  /*
     Takes in a row of G and returns indexes of coeff to keep in the row sparse version
     energy:: how much energy - in l2 sens - to keep after hard-thresholding
     */
 
  if (energy >= 1) return 0;
 
  const Sparse<t_real> row_abs = row.cwiseAbs();
 
  t_real abs_row_max = 0;
  for (t_uint i = 0; i < row_abs.nonZeros(); i++) {
    if (*(row_abs.valuePtr() + i) > abs_row_max) abs_row_max = *(row_abs.valuePtr() + i);
  }
  const t_real row_total_energy = row_abs.cwiseProduct(row_abs).sum();
  const t_real target_threshold_energy_sum = row_total_energy * energy;
  const t_int niters = 200;
  t_real lower = 0;
  t_real upper = abs_row_max;
  t_real current_threshold = upper * 0.5;
  t_real old_threshold = upper;
  bool converged = false;
  for (t_int i = 0; i < niters; i++) {
    t_real threshold_energy_sum = 0;
    for (t_int j = 0; j < row_abs.nonZeros(); j++) {
      if (*(row_abs.valuePtr() + j) > current_threshold)
        threshold_energy_sum += *(row_abs.valuePtr() + j) * *(row_abs.valuePtr() + j);
    }
    if ((std::abs(current_threshold - old_threshold) / std::abs(old_threshold) < 1e-3) and
        ((threshold_energy_sum >= target_threshold_energy_sum) and
         (((threshold_energy_sum / target_threshold_energy_sum) - 1) < 0.01))) {
      converged = true;
      break;
    }
    if (threshold_energy_sum < target_threshold_energy_sum)
      upper = current_threshold;
    else
      lower = current_threshold;
    old_threshold = current_threshold;
    current_threshold = (lower + upper) * 0.5;
  }
  if (!converged) current_threshold = lower;
  return current_threshold;
}

Referenced by wprojection_matrix().

sparsity_im()

t_real purify::wproj_utilities::sparsity_im

(

const Image< t_complex > &

Cmat

)

return faction of non zero values from matrix

sparsity_sp()

t_real purify::wproj_utilities::sparsity_sp

(

const Sparse< t_complex > &

Gmat

)

return fraction of non zero values from sparse matrix

upsample_ratio_sim()

t_real purify::wproj_utilities::upsample_ratio_sim	(	const utilities::vis_params &	uv_vis,
		const t_real &	L,
		const t_real &	M,
		const t_int &	x_size,
		const t_int &	y_size,
		const t_int &	multipleOf
	)

Calculate upsample ratio from bandwidth (only needed for simulations)

w_rows()

std::vector<t_uint> purify::wproj_utilities::w_rows

(

const Sparse< t_complex > &

w_degrider

)

Definition at line 22 of file wproj_utilities.cc.

                                                              {
  std::set<t_uint> w_rows_set;
  for (t_int k = 0; k < w_degrider.outerSize(); ++k) {
    t_uint index = 0;
    for (Sparse<t_complex>::InnerIterator it(w_degrider, k); it; ++it) {
      index++;
      if (index > 0) {
        w_rows_set.insert(it.col());
      }
    }
  }
  std::vector<t_uint> w_rows(w_rows_set.begin(), w_rows_set.end());
  std::sort(w_rows.begin(), w_rows.end());
  return w_rows;
}

wprojection_matrix()

Sparse< t_complex > purify::wproj_utilities::wprojection_matrix	(	const Sparse< t_complex > &	G,
		const t_uint &	x_size,
		const t_uint &	y_size,
		const Vector< t_real > &	w_components,
		const t_real &	cell_x,
		const t_real &	cell_y,
		const t_real &	energy_fraction_chirp,
		const t_real &	energy_fraction_wproj,
		const expansions::series	series,
		const t_uint	order,
		const t_real &	interpolation_error
	)

Produce Gridding matrix convovled with chirp matrix for wprojection.

Definition at line 132 of file wproj_utilities.cc.

                                                                        {
  const t_uint Npix = x_size * y_size;
  const t_uint Nvis = w_components.size();
 
  PURIFY_HIGH_LOG("Spread of w components {} ",
                  std::sqrt(std::pow(w_components.norm(), 2) / Nvis -
                            std::pow(w_components.array().mean(), 2)));
  PURIFY_HIGH_LOG("Mean of w components {} ", w_components.array().mean());
  PURIFY_HIGH_LOG("Hard-thresholding of the chirp kernels: energy [{}] ", energy_fraction_chirp);
  PURIFY_HIGH_LOG("Hard-thresholding of the rows of G: energy [{}]  ", energy_fraction_wproj);
 
  const auto ft_plan = operators::fftw_plan::measure;
  const auto fftop_ = operators::init_FFT_2d<Vector<t_complex>>(y_size, x_size, 1., ft_plan);
  Sparse<t_complex> GW(Nvis, Npix);
 
  t_uint counts = 0;
  t_uint total_non_zero = 0;
  GW.reserve(G.nonZeros() * 1e3);
#pragma omp parallel for
  for (t_int m = 0; m < G.rows(); m++) {
    const Sparse<t_complex> chirp =
        create_chirp_row(w_components(m), cell_x, cell_y, x_size, y_size, energy_fraction_chirp,
                         std::get<0>(fftop_));
 
    Sparse<t_complex> kernel = row_wise_sparse_convolution(G.row(m), chirp, x_size, y_size);
    const t_real thres = sparsify_row_thres(kernel, energy_fraction_wproj);
    kernel.prune([&](const t_uint &i, const t_uint &j, const t_complex &value) {
      return std::abs(value) > thres;
    });
    assert(kernel.cols() > 0);
    assert(kernel.rows() == 1);
#pragma omp critical
    GW.row(m) = kernel;
#pragma omp critical
    counts++;
#pragma omp critical
    total_non_zero += kernel.nonZeros();
#pragma omp critical
    // if(counts % 100 == 0) {
    PURIFY_DEBUG("Row {} of {}", counts, GW.rows());
    PURIFY_DEBUG("With {} entries non zero, which is {} entries per a row.", total_non_zero,
                 static_cast<t_real>(total_non_zero) / counts);
    //}
  }
  assert(GW.nonZeros() > 0);
  PURIFY_DEBUG("\nBuilding the rows of GW.. DONE!\n");
  PURIFY_DEBUG("DONE - With {} entries non zero, which is {} entries per a row.", GW.nonZeros(),
               static_cast<t_real>(GW.nonZeros()) / GW.rows());
  GW.makeCompressed();
  return GW;
}

References create_chirp_row(), purify::operators::measure, PURIFY_DEBUG, PURIFY_HIGH_LOG, row_wise_sparse_convolution(), and sparsify_row_thres().