lsst::gauss2d::detail Namespace Reference

Classes
class	GradientsExtra
struct	Terms
class	TermsGradient
struct	TermsMoments
class	TermsPixel
	Storage for terms common to Gaussians for a single pixel. More...
struct	ValuesGauss

Typedefs
typedef std::vector< double >	vecd
typedef std::unique_ptr< vecd >	vecdptr
typedef std::array< double, N_PARAMS_GAUSS2D >	Weights
typedef std::vector< TermsPixel >	TermsPixelVec
typedef std::vector< TermsGradient >	TermsGradientVec
using	type_name_prober = void

Functions
TermsMoments	gaussian_pixel_x_xx (const double cen_x, const double x_min, const double bin_x, const unsigned int dim_x, const double xx_weight, const double xy_weight)
Terms	terms_from_covar (double weight, const Ellipse &ell)
template<class Data>
void	gaussian_pixel (Data &image, const double weight, const double cen_x, const double cen_y, const double xx_weight, const double yy_weight, const double xy_weight)
void	gaussian_pixel_get_jacobian (ValuesGauss &out, const double m, const double m_unweight, const double xmc_norm, const double ymc_weighted, const double xmc, const double ymc, const double norms_yy, const double xmc_t_xy_factor, const double sig_x_inv, const double sig_y_inv, const double xy_norm, const double xx_norm, const double yy_weighted, const double rho_div_one_m_rhosq, const double norms_xy_div_rho, const double dsig_x_conv_src=1, const double dsig_y_conv_src=1, const double drho_c_dsig_x_src=0, const double drho_c_dsig_y_src=0, const double drho_c_drho_s=1)
void	gaussian_pixel_get_jacobian_from_terms (ValuesGauss &out, const size_t dim1, const TermsPixel &terms_pixel, const TermsGradient &terms_grad, double m, double m_unweighted, double xy_norm)
template<typename t>
void	gaussian_pixel_add_values (t &cen_x, t &cen_y, t &L, t &sig_x, t &sig_y, t &rho, const ValuesGauss &values, const t weight_cen_x=1, const t weight_cen_y=1, const t weight_L=1, const t weight_sig_x=1, const t weight_sig_y=1, const t weight_rho=1)
template<typename t, class Data, class Indices, bool do_extra>
void	gaussians_pixel_add_like_grad (Image< t, Data > *output, const Image< idx_type, Indices > &grad_param_map, const Image< t, Data > &grad_param_factor, const size_t N_GAUSS, const std::vector< Weights > &gaussweights, double &chi_pix, double &loglike, const double model, double data, double sigma_inv, unsigned int dim1, unsigned int dim2, const TermsPixelVec &terms_pixel, const TermsGradientVec &terms_grad, ValuesGauss &gradients, const std::shared_ptr< const Image< idx_type, Indices > > extra_param_map, const std::shared_ptr< const Image< t, Data > > extra_param_factor)
template<typename T, class Data, class Indices, GradientType gradient_type, bool do_extra>
T	gaussian_pixel_add_all (size_t g, size_t j, size_t i, double weight, double sigma_inv, const TermsPixelVec &terms_pixel_vec, ImageArray< T, Data > &output_jac_ref, const Image< idx_type, Indices > &grad_param_map, const Image< T, Data > &grad_param_factor, std::vector< Weights > &gradweights, const TermsGradientVec &terms_grad_vec, ValuesGauss &gradients, GradientsExtra< T, Data, Indices > *grad_extra)
template<typename Indices>
const std::shared_ptr< const Indices >	_param_map_default (size_t n_gaussians, size_t n_params=N_PARAMS_GAUSS2D, size_t increment=1)
template<typename Data>
const std::shared_ptr< const Data >	_param_factor_default (size_t n_gaussians, size_t n_params=N_PARAMS_GAUSS2D, double value=1.)
template<typename T>
constexpr std::string_view	wrapped_type_name ()
constexpr std::size_t	wrapped_type_name_prefix_length ()
constexpr std::size_t	wrapped_type_name_suffix_length ()

Variables
constexpr std::string_view	NAMESPACE_SEPARATOR = "::"

Typedef Documentation

TermsGradientVec

typedef std::vector<TermsGradient> lsst::gauss2d::detail::TermsGradientVec

Definition at line 327 of file evaluate.h.

TermsPixelVec

typedef std::vector<TermsPixel> lsst::gauss2d::detail::TermsPixelVec

Definition at line 229 of file evaluate.h.

type_name_prober

using lsst::gauss2d::detail::type_name_prober = void

Definition at line 47 of file type_name.h.

vecd

typedef std::vector<double> lsst::gauss2d::detail::vecd

Definition at line 65 of file evaluate.h.

vecdptr

typedef std::unique_ptr<vecd> lsst::gauss2d::detail::vecdptr

Definition at line 66 of file evaluate.h.

Weights

typedef std::array<double, N_PARAMS_GAUSS2D> lsst::gauss2d::detail::Weights

Definition at line 179 of file evaluate.h.

Function Documentation

_param_factor_default()

template<typename Data>

const std::shared_ptr< const Data > lsst::gauss2d::detail::_param_factor_default	(	size_t	n_gaussians,
		size_t	n_params = N_PARAMS_GAUSS2D,
		double	value = 1. )

Definition at line 644 of file evaluate.h.

                                                                           {
    auto param_factor = std::make_shared<Data>(n_gaussians, n_params);
    for (size_t g = 0; g < n_gaussians; ++g) {
        for (size_t p = 0; p < n_params; ++p) {
            param_factor->set_value(g, p, value);
        }
    }
    return std::const_pointer_cast<const Data>(param_factor);
}

_param_map_default()

template<typename Indices>

const std::shared_ptr< const Indices > lsst::gauss2d::detail::_param_map_default	(	size_t	n_gaussians,
		size_t	n_params = N_PARAMS_GAUSS2D,
		size_t	increment = 1 )

Definition at line 629 of file evaluate.h.

                                                                              {
    auto param_map = std::make_shared<Indices>(n_gaussians, n_params);
    size_t index = 0;
    for (size_t g = 0; g < n_gaussians; ++g) {
        for (size_t p = 0; p < n_params; ++p) {
            param_map->set_value(g, p, index);
            index += increment;
        }
    }
    return std::const_pointer_cast<const Indices>(param_map);
}

gaussian_pixel()

template<class Data>

void lsst::gauss2d::detail::gaussian_pixel ( Data & image, const double weight, const double cen_x, const double cen_y, const double xx_weight, const double yy_weight, const double xy_weight )

inline

Definition at line 407 of file evaluate.h.

                                                                                                   {
    const unsigned int dim_x = image.get_n_cols();
    const unsigned int dim_y = image.get_n_rows();
    const auto coordsys = image.get_coordsys();
    const double x_min = coordsys.x_min;
    const double y_min = coordsys.y_min;
    const double bin_x = coordsys.get_dx1() / image.get_n_cols();
    const double bin_y = coordsys.get_dy2() / image.get_n_rows();
 
    const double weight_pix = weight * bin_x * bin_y;
 
    const auto moment_terms_y = gaussian_pixel_x_xx(cen_y, y_min, bin_y, dim_y, yy_weight, xy_weight);
    const vecd& y_weighted = *(moment_terms_y.x_norm);
    const vecd& yy_weighted = *(moment_terms_y.xx);
 
    double x = x_min - cen_x + bin_x / 2.;
    // TODO: Consider a version with cached xy, although it doubles memory usage
    for (unsigned int i = 0; i < dim_x; i++) {
        const double xx_weighted = x * x * xx_weight;
        for (unsigned int j = 0; j < dim_y; j++) {
            image.set_value_unchecked(j, i,
                                      weight_pix * exp(-(xx_weighted + yy_weighted[j] - x * y_weighted[j])));
        }
        x += bin_x;
    }
}

gaussian_pixel_add_all()

template<typename T, class Data, class Indices, GradientType gradient_type, bool do_extra>

T lsst::gauss2d::detail::gaussian_pixel_add_all ( size_t g, size_t j, size_t i, double weight, double sigma_inv, const TermsPixelVec & terms_pixel_vec, ImageArray< T, Data > & output_jac_ref, const Image< idx_type, Indices > & grad_param_map, const Image< T, Data > & grad_param_factor, std::vector< Weights > & gradweights, const TermsGradientVec & terms_grad_vec, ValuesGauss & gradients, GradientsExtra< T, Data, Indices > * grad_extra )

inline

Definition at line 592 of file evaluate.h.

                                                                              {
    const TermsPixel& terms_pixel = terms_pixel_vec[g];
    const double xy_norm = terms_pixel.xmc * (*terms_pixel.ymc_weighted)[j];
    const double value_unweight
            = terms_pixel.weight * exp(-(terms_pixel.xmc_sq_norm + (*terms_pixel.yy_weighted)[j] - xy_norm));
    const double value = weight * value_unweight;
    // Computes dmodel/dx for x in [cen_x, cen_y, flux, sigma_x, sigma_y, rho]
    if constexpr (gradient_type == GradientType::loglike) {
        Weights& output = gradweights[g];
        output[0] = value;
        output[1] = value_unweight;
        output[2] = xy_norm;
    } else if constexpr (gradient_type == GradientType::jacobian) {
        gaussian_pixel_get_jacobian_from_terms(gradients, j, terms_pixel, terms_grad_vec[g],
                                               sigma_inv * value, sigma_inv * value_unweight, xy_norm);
        gaussian_pixel_add_values<T>(
                output_jac_ref[grad_param_map.get_value_unchecked(g, 0)]._get_value_unchecked(j, i),
                output_jac_ref[grad_param_map.get_value_unchecked(g, 1)]._get_value_unchecked(j, i),
                output_jac_ref[grad_param_map.get_value_unchecked(g, 2)]._get_value_unchecked(j, i),
                output_jac_ref[grad_param_map.get_value_unchecked(g, 3)]._get_value_unchecked(j, i),
                output_jac_ref[grad_param_map.get_value_unchecked(g, 4)]._get_value_unchecked(j, i),
                output_jac_ref[grad_param_map.get_value_unchecked(g, 5)]._get_value_unchecked(j, i),
                gradients, grad_param_factor.get_value_unchecked(g, 0),
                grad_param_factor.get_value_unchecked(g, 1), grad_param_factor.get_value_unchecked(g, 2),
                grad_param_factor.get_value_unchecked(g, 3), grad_param_factor.get_value_unchecked(g, 4),
                grad_param_factor.get_value_unchecked(g, 5));
        if constexpr (do_extra) grad_extra->add(g, j, i, gradients);
    }
    return value;
}

gaussian_pixel_add_values()

template<typename t>

void lsst::gauss2d::detail::gaussian_pixel_add_values ( t & cen_x, t & cen_y, t & L, t & sig_x, t & sig_y, t & rho, const ValuesGauss & values, const t weight_cen_x = 1, const t weight_cen_y = 1, const t weight_L = 1, const t weight_sig_x = 1, const t weight_sig_y = 1, const t weight_rho = 1 )

inline

Definition at line 533 of file evaluate.h.

                                                              {
    cen_x += weight_cen_x * values.cen_x;
    cen_y += weight_cen_y * values.cen_y;
    L += weight_L * values.L;
    sig_x += weight_sig_x * values.sigma_x;
    sig_y += weight_sig_y * values.sigma_y;
    rho += weight_rho * values.rho;
}

gaussian_pixel_get_jacobian()

void lsst::gauss2d::detail::gaussian_pixel_get_jacobian ( ValuesGauss & out, const double m, const double m_unweight, const double xmc_norm, const double ymc_weighted, const double xmc, const double ymc, const double norms_yy, const double xmc_t_xy_factor, const double sig_x_inv, const double sig_y_inv, const double xy_norm, const double xx_norm, const double yy_weighted, const double rho_div_one_m_rhosq, const double norms_xy_div_rho, const double dsig_x_conv_src = 1, const double dsig_y_conv_src = 1, const double drho_c_dsig_x_src = 0, const double drho_c_dsig_y_src = 0, const double drho_c_drho_s = 1 )

inline

Definition at line 495 of file evaluate.h.

                                                                                                            {
    out.cen_x = m * (2 * xmc_norm - ymc_weighted);
    out.cen_y = m * (2 * ymc * norms_yy - xmc_t_xy_factor);
    out.L = m_unweight;
    const double one_p_xy = 1. + xy_norm;
    // df/dsigma_src = df/d_sigma_conv * dsigma_conv/dsigma_src
    const double dfdrho_c = m
                            * (rho_div_one_m_rhosq * (1 - 2 * (xx_norm + yy_weighted - xy_norm))
                               + xmc * ymc * norms_xy_div_rho);
    out.sigma_x = dfdrho_c * drho_c_dsig_x_src + dsig_x_conv_src * (m * sig_x_inv * (2 * xx_norm - one_p_xy));
    out.sigma_y
            = dfdrho_c * drho_c_dsig_y_src + dsig_y_conv_src * (m * sig_y_inv * (2 * yy_weighted - one_p_xy));
    // drho_conv/drho_src = sigmaxy_src/sigmaxy_conv
    out.rho = dfdrho_c * drho_c_drho_s;
}

gaussian_pixel_get_jacobian_from_terms()

void lsst::gauss2d::detail::gaussian_pixel_get_jacobian_from_terms ( ValuesGauss & out, const size_t dim1, const TermsPixel & terms_pixel, const TermsGradient & terms_grad, double m, double m_unweighted, double xy_norm )

inline

Definition at line 519 of file evaluate.h.

                                                                                        {
    gaussian_pixel_get_jacobian(
            out, m, m_unweighted * terms_pixel.weight_kernel, terms_pixel.xmc_weighted,
            (*terms_pixel.ymc_weighted)[dim1], terms_pixel.xmc, (*terms_grad.ymc)[dim1], terms_grad.yy_weight,
            terms_grad.xmc_t_xy_weight, terms_grad.sig_x_inv, terms_grad.sig_y_inv, xy_norm,
            terms_pixel.xmc_sq_norm, (*terms_pixel.yy_weighted)[dim1], terms_grad.rho_factor,
            terms_grad.rho_xy_factor, terms_grad.sig_x_src_div_conv, terms_grad.sig_y_src_div_conv,
            terms_grad.drho_c_dsig_x_src, terms_grad.drho_c_dsig_y_src, terms_grad.drho_c_drho_s);
}

gaussian_pixel_x_xx()

TermsMoments lsst::gauss2d::detail::gaussian_pixel_x_xx ( const double cen_x, const double x_min, const double bin_x, const unsigned int dim_x, const double xx_weight, const double xy_weight )

inline

Definition at line 83 of file evaluate.h.

                                                                {
    const double x_init = x_min - cen_x + bin_x / 2.;
    vecdptr x = std::make_unique<vecd>(dim_x);
    vecdptr xx = std::make_unique<vecd>(dim_x);
    vecdptr x_norm = std::make_unique<vecd>(dim_x);
    for (unsigned int i = 0; i < dim_x; i++) {
        double dist = x_init + i * bin_x;
        (*x)[i] = dist;
        (*xx)[i] = dist * dist * xx_weight;
        (*x_norm)[i] = dist * xy_weight;
    }
    return TermsMoments({std::move(x), std::move(x_norm), std::move(xx)});
}

gaussians_pixel_add_like_grad()

template<typename t, class Data, class Indices, bool do_extra>

void lsst::gauss2d::detail::gaussians_pixel_add_like_grad ( Image< t, Data > * output, const Image< idx_type, Indices > & grad_param_map, const Image< t, Data > & grad_param_factor, const size_t N_GAUSS, const std::vector< Weights > & gaussweights, double & chi_pix, double & loglike, const double model, double data, double sigma_inv, unsigned int dim1, unsigned int dim2, const TermsPixelVec & terms_pixel, const TermsGradientVec & terms_grad, ValuesGauss & gradients, const std::shared_ptr< const Image< idx_type, Indices > > extra_param_map, const std::shared_ptr< const Image< t, Data > > extra_param_factor )

inline

Definition at line 548 of file evaluate.h.

                                                                      {
    double diff = data - model;
    chi_pix = sigma_inv * diff;
    double diffvar = sigma_inv * chi_pix;
    loglike -= diff * diffvar / 2.;
    for (size_t g = 0; g < N_GAUSS; ++g) {
        const Weights& weights = gaussweights[g];
        /*
        * Derivation:
        *
            ll = sum(-(data-model)^2*varinv/2)
            dll/dx = --2*dmodel/dx*(data-model)*varinv/2
            dmodelsum/dx = d(.. + model[g])/dx = dmodel[g]/dx
            dll/dx = dmodel[g]/dx*diffvar
        */
        gaussian_pixel_get_jacobian_from_terms(gradients, dim1, terms_pixel[g], terms_grad[g],
                                               weights[0] * diffvar, weights[1] * diffvar, weights[2]);
        gaussian_pixel_add_values<t>(
                output->_get_value_unchecked(0, grad_param_map.get_value_unchecked(g, 0)),
                output->_get_value_unchecked(0, grad_param_map.get_value_unchecked(g, 1)),
                output->_get_value_unchecked(0, grad_param_map.get_value_unchecked(g, 2)),
                output->_get_value_unchecked(0, grad_param_map.get_value_unchecked(g, 3)),
                output->_get_value_unchecked(0, grad_param_map.get_value_unchecked(g, 4)),
                output->_get_value_unchecked(0, grad_param_map.get_value_unchecked(g, 5)), gradients,
                grad_param_factor.get_value_unchecked(g, 0), grad_param_factor.get_value_unchecked(g, 1),
                grad_param_factor.get_value_unchecked(g, 2), grad_param_factor.get_value_unchecked(g, 3),
                grad_param_factor.get_value_unchecked(g, 4), grad_param_factor.get_value_unchecked(g, 5));
        if constexpr (do_extra) {
            double value = gradients.L * extra_param_factor->get_value_unchecked(g, 0)
                           + gradients.sigma_x * extra_param_factor->get_value_unchecked(g, 1)
                           + gradients.sigma_y * extra_param_factor->get_value_unchecked(g, 2);
            output->_get_value_unchecked(0, extra_param_map->get_value_unchecked(g, 1)) += value;
        }
    }
}

terms_from_covar()

Terms lsst::gauss2d::detail::terms_from_covar	(	double	weight,
		const Ellipse &	ell )

Definition at line 30 of file evaluate.cc.

                                                                {
    double sig_x = ell.get_sigma_x();
    double sig_y = ell.get_sigma_y();
    double rho = ell.get_rho();
 
    const double norm_exp = 1. / (1 - rho * rho);
    Terms rval = {.weight = weight / (2. * M_PI * sig_x * sig_y) * sqrt(norm_exp),
                  .xx = norm_exp / (2. * sig_x * sig_x),
                  .yy = norm_exp / (2. * sig_y * sig_y),
                  .xy = rho * norm_exp / (sig_x * sig_y)};
    return rval;
}

wrapped_type_name()

template<typename T>

std::string_view lsst::gauss2d::detail::wrapped_type_name ( )

constexpr

Definition at line 50 of file type_name.h.

                                             {
#ifdef __clang__
    return __PRETTY_FUNCTION__;
#elif defined(__GNUC__)
    return __PRETTY_FUNCTION__;
#elif defined(_MSC_VER)
    return __FUNCSIG__;
#else
#error "Unsupported compiler"
#endif
}

wrapped_type_name_prefix_length()

std::size_t lsst::gauss2d::detail::wrapped_type_name_prefix_length ( )

constexpr

Definition at line 62 of file type_name.h.

                                                      {
    return wrapped_type_name<type_name_prober>().find(type_name<type_name_prober>());
}

wrapped_type_name_suffix_length()

std::size_t lsst::gauss2d::detail::wrapped_type_name_suffix_length ( )

constexpr

Definition at line 66 of file type_name.h.

                                                      {
    return wrapped_type_name<type_name_prober>().length() - wrapped_type_name_prefix_length()
           - type_name<type_name_prober>().length();
}

Variable Documentation

NAMESPACE_SEPARATOR

std::string_view lsst::gauss2d::detail::NAMESPACE_SEPARATOR = "::"

constexpr

Definition at line 71 of file type_name.h.