lsst::gauss2d::fit::Model< T, Image, Indices, Mask > Class Template Reference

A collection of Sources comprising a ParametricModel of Data. More...

#include <model.h>

Inheritance diagram for lsst::gauss2d::fit::Model< T, Image, Indices, Mask >:

Public Types
typedef GaussianEvaluator< T, Image, Indices >	Evaluator
typedef std::map< std::reference_wrapper< const Channel >, std::vector< std::unique_ptr< const Gaussians > > >	GaussiansMap
typedef Data< T, Image, Mask >	ModelData
typedef ModelData::Observation	Observation

Public Member Functions
	Model (std::shared_ptr< const ModelData > data, PsfModels &psfmodels, Sources &sources, Priors &priors)
	Construct a Model instance from Data, PsfModels, Sources and Priors.
void	add_extra_param_map (const Channel &channel, ExtraParamMap &map_extra, const GradParamMap &map_grad, ParameterMap &offsets) const override
	Add extra Parameter indices to a map.
void	add_extra_param_factors (const Channel &channel, ExtraParamFactors &factors) const override
	Add extra Parameter gradient factors to an existing vector.
void	add_grad_param_map (const Channel &channel, GradParamMap &map, ParameterMap &offsets) const override
	Add Parameter gradient indices to an existing map.
void	add_grad_param_factors (const Channel &channel, GradParamFactors &factors) const override
	Add Parameter gradient factors to an existing map.
std::vector< double >	compute_loglike_grad (bool include_prior=true, bool print=false, bool verify=false, double findiff_frac=1e-5, double findiff_add=1e-5, double rtol=1e-3, double atol=1e-8)
	Compute the gradient (partial first derivative) of the log-likelihood for each free parameter.
std::unique_ptr< Image >	compute_hessian (bool transformed=false, bool include_prior=true, std::optional< HessianOptions > options=std::nullopt, bool print=false)
	Compute the Hessian matrix (second order partial derivatives) of the log likehood.
std::vector< double >	evaluate (bool print=false, bool normalize_loglike=false)
	Evaluate the model for every Observation in _data.
double	evaluate_observation (size_t idx)
	Evaluate a single observation with the given index in _data.
std::shared_ptr< const ModelData >	get_data () const
	Return _data.
EvaluatorMode	get_mode () const
std::unique_ptr< const lsst::gauss2d::Gaussians >	get_gaussians (const Channel &channel) const override
	Return the vector of Gaussian sub-components controlled by this model.
std::vector< double >	get_loglike_const_terms ()
	Get the constant (variance-dependent) terms of the log likelihood for each observation.
size_t	get_n_gaussians (const Channel &channel) const override
	Return the number of Gaussian sub-components controlled by this model.
std::vector< std::shared_ptr< Image > >	get_outputs () const
	Return _outputs (output Image instances for each Observation in _data)
std::vector< std::pair< ParamBaseCRef, size_t > >	get_offsets_parameters () const
ParamRefs &	get_parameters (ParamRefs &params, ParamFilter *filter=nullptr) const override
	Add Parameter refs matching the filter to a vector, in order.
ParamCRefs &	get_parameters_const (ParamCRefs &params, ParamFilter *filter=nullptr) const override
	Same as get_parameters(), but for const refs.
ParamRefs &	get_parameters_observation (ParamRefs &params, size_t idx, ParamFilter *filter=nullptr) const
	Same as get_parameters(), but for a single Observation with index idx in _data.
ParamCRefs &	get_parameters_observation_const (ParamCRefs &params, size_t idx, ParamFilter *filter=nullptr) const
	Same as get_parameters_const(), but for a single Observation with index idx in _data.
Priors	get_priors () const
	Return _priors, the list of Prior instances.
PsfModels	get_psfmodels () const
	Return _psfmodels, the list of PsfModel instances for each Observation in _data.
Sources	get_sources () const
	Return _sources, the list of Source instances for each Observation in _data.
void	set_extra_param_factors (const Channel &channel, ExtraParamFactors &factors, size_t index) const override
	Set extra Parameter gradient factors in an existing map.
void	set_grad_param_factors (const Channel &channel, GradParamFactors &factors, size_t index) const override
	Set Parameter gradient factors in an existing map.
void	setup_evaluators (EvaluatorMode mode=EvaluatorMode::image, std::vector< std::vector< std::shared_ptr< Image > > > outputs={}, std::vector< std::shared_ptr< Image > > residuals={}, std::vector< std::shared_ptr< Image > > outputs_prior={}, std::shared_ptr< Image > residuals_prior=nullptr, bool force=false, bool print=false)
	Setup Evaluator instances for every Observation in _data using the given EvaluatorMode.
size_t	size () const
	Get the size of this->_data.
std::string	repr (bool name_keywords=false, std::string_view namespace_separator=Object::CC_NAMESPACE_SEPARATOR) const override
	Return a full, callable string representation of this.
std::string	str () const override
	Return a brief, human-readable string representation of this.
std::vector< std::string >	verify_jacobian (double findiff_frac=1e-5, double findiff_add=1e-5, double rtol=1e-3, double atol=1e-3, double max_ll_diff=0)
	Verify that the Jacobian is correct by comparing to finite differences.
ParamRefs	get_parameters_new (ParamFilter *filter=nullptr) const
	Same as get_parameters(), but returning a new vector.
ParamCRefs	get_parameters_const_new (ParamFilter *filter=nullptr) const
	Same as get_parameters_const(), but returning a new vector.

Static Public Member Functions
static std::string_view	null_str (const std::string_view &namespace_separator)

Static Public Attributes
static constexpr std::string_view	CC_NAMESPACE_SEPARATOR = "::"
	The C++ namespace separator.
static constexpr std::string_view	NULL_STR_GENERAL = "None"
static constexpr std::string_view	PY_NAMESPACE_SEPARATOR = "."

Detailed Description

template<typename T, typename Image, typename Indices, typename Mask>
class lsst::gauss2d::fit::Model< T, Image, Indices, Mask >

A collection of Sources comprising a ParametricModel of Data.

Template Parameters

T	The type of the Image (usually float or double).
Image	The class of image used in Data.
Indices	The class of index map used in Data.
Mask	The class of mask used in Data.

Definition at line 78 of file model.h.

Member Typedef Documentation

Evaluator

template<typename T, typename Image, typename Indices, typename Mask>

typedef GaussianEvaluator<T, Image, Indices> lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::Evaluator

Definition at line 80 of file model.h.

GaussiansMap

template<typename T, typename Image, typename Indices, typename Mask>

typedef std::map<std::reference_wrapper<const Channel>, std::vector<std::unique_ptr<const Gaussians> > > lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::GaussiansMap

Definition at line 82 of file model.h.

ModelData

template<typename T, typename Image, typename Indices, typename Mask>

typedef Data<T, Image, Mask> lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::ModelData

Definition at line 83 of file model.h.

Observation

template<typename T, typename Image, typename Indices, typename Mask>

typedef ModelData::Observation lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::Observation

Definition at line 84 of file model.h.

Constructor & Destructor Documentation

Model()

template<typename T, typename Image, typename Indices, typename Mask>

lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::Model ( std::shared_ptr< const ModelData > data, PsfModels & psfmodels, Sources & sources, Priors & priors )

inlineexplicit

Construct a Model instance from Data, PsfModels, Sources and Priors.

Parameters

data	The data to model.
psfmodels	A vector of PSF models, ordered to match each Observation in data.
sources	A vector of Source models.
priors	A vector of Prior likelihoods.

Definition at line 94 of file model.h.

            : _data(std::move(data)),
              _size(_data == nullptr ? 0 : _data->size()),
              _size_priors(priors.size()) {
        if (_data == nullptr) throw std::invalid_argument("Model data can't be null");
        size_t size = this->size();
        if (psfmodels.size() != size) {
            throw std::invalid_argument("Model psfmodels.size()=" + std::to_string(psfmodels.size())
                                        + "!= data.size()=" + std::to_string(size));
        }
 
        _outputs.resize(size);
        _psfmodels.reserve(size);
        size_t i = 0;
        for (auto& psfmodel : psfmodels) {
            if (psfmodel == nullptr)
                throw std::invalid_argument("Model psfmodels[" + std::to_string(i) + "] can't be null");
            _psfcomps.push_back(psfmodel->get_gaussians());
            _psfmodels.push_back(psfmodel);
            i++;
        }
 
        _priors.reserve(_size_priors);
        i = 0;
        for (auto& prior : priors) {
            if (prior == nullptr)
                throw std::invalid_argument("Model priors[" + std::to_string(i) + "] can't be null");
            _priors.push_back(prior);
            i++;
        }
 
        _sources.reserve(sources.size());
        i = 0;
        for (auto& source : sources) {
            if (source == nullptr)
                throw std::invalid_argument("Model Sources[" + std::to_string(i) + "] can't be null");
            _sources.push_back(source);
            i++;
        }
 
        _gaussians_srcs = this->_get_gaussians_srcs();
        _factors_extra_in.resize(_size);
        _factors_grad_in.resize(_size);
        _map_extra_in.resize(_size);
        _map_grad_in.resize(_size);
    }

Member Function Documentation

add_extra_param_factors()

template<typename T, typename Image, typename Indices, typename Mask>

void lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::add_extra_param_factors ( const Channel & channel, ExtraParamFactors & factors ) const

inlineoverridevirtual

Add extra Parameter gradient factors to an existing vector.

Parameters

channel	The Channel to add factors for.
factors	The ExtraParamFactors to add to.

Implements lsst::gauss2d::fit::ParametricModel.

Definition at line 681 of file model.h.

                                                                                                    {
        for (auto& source : _sources) source->add_extra_param_factors(channel, factors);
    }

add_extra_param_map()

template<typename T, typename Image, typename Indices, typename Mask>

void lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::add_extra_param_map ( const Channel & channel, ExtraParamMap & map_extra, const GradParamMap & map_grad, ParameterMap & offsets ) const

inlineoverridevirtual

Add extra Parameter indices to a map.

Parameters

channel	The Channel to add indices for.
map_extra	The ExtraParamMap to add to.
map_grad	The completed GradParamMap.
offsets	A map of index offsets for Parameters that have already been added.

Implements lsst::gauss2d::fit::ParametricModel.

Definition at line 677 of file model.h.

                                                                   {
        for (auto& source : _sources) source->add_extra_param_map(channel, map_extra, map_grad, offsets);
    }

add_grad_param_factors()

template<typename T, typename Image, typename Indices, typename Mask>

void lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::add_grad_param_factors ( const Channel & channel, GradParamFactors & factors ) const

inlineoverridevirtual

Add Parameter gradient factors to an existing map.

Parameters

channel	The Channel to add factors for.
factors	The GradParamFactors to add to.

Implements lsst::gauss2d::fit::ParametricModel.

Definition at line 687 of file model.h.

                                                                                                  {
        for (auto& source : _sources) source->add_grad_param_factors(channel, factors);
    }

add_grad_param_map()

template<typename T, typename Image, typename Indices, typename Mask>

void lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::add_grad_param_map ( const Channel & channel, GradParamMap & map, ParameterMap & offsets ) const

inlineoverridevirtual

Add Parameter gradient indices to an existing map.

Parameters

channel	The Channel to add indices for.
map	The map to add to.
offsets	A map of index offsets for Parameters that have already been added.

Implements lsst::gauss2d::fit::ParametricModel.

Definition at line 684 of file model.h.

                                                                                                             {
        for (auto& source : _sources) source->add_grad_param_map(channel, map, offsets);
    }

compute_hessian()

template<typename T, typename Image, typename Indices, typename Mask>

std::unique_ptr< Image > lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::compute_hessian ( bool transformed = false, bool include_prior = true, std::optional< HessianOptions > options = std::nullopt, bool print = false )

inline

Compute the Hessian matrix (second order partial derivatives) of the log likehood.

Parameters

transformed	Whether the matrix should be computed for transformed parameters or not If not, parameter transforms are temporarily removed.
include_prior	Whether to include the prior likelihood(s) in the Hessian.
options	Options for computing the Hessian via finite differencing of loglikelihood gradients. If null, the Hessian is estimated as J^T J (where J is the Jacobian).

Returns

The Hessian matrix for all free parameters.

Definition at line 825 of file model.h.

                                                               {
        this->setup_evaluators(EvaluatorMode::loglike_grad, {}, {}, {}, nullptr, true, print);
 
        auto filter_free = g2f::ParamFilter{false, true, true, true};
        auto params_free = this->get_parameters_new(&filter_free);
        params_free = nonconsecutive_unique<ParamBaseRef>(params_free);
        const size_t n_params_free = params_free.size();
 
        std::vector<std::shared_ptr<const parameters::Transform<double>>> transforms = {};
        transforms.reserve(n_params_free);
        if (!transformed) {
            for (auto& paramref : params_free) {
                auto& param = paramref.get();
                transforms.emplace_back(param.get_transform_ptr());
                double value = param.get_value();
                param.set_transform(nullptr);
                if (param.get_value() != value) {
                    throw std::logic_error("Param " + param.str() + " changed value from "
                                           + std::to_string(value) + " after dropping transform");
                }
            }
        }
 
        auto hessian = std::make_unique<Image>(n_params_free, n_params_free);
        const auto n_obs = this->size();
        std::vector<size_t> param_idx(n_params_free);
 
        size_t n_integral = 0;
        size_t n_frac = 0;
 
        for (size_t idx_param = 0; idx_param < n_params_free; idx_param++) {
            const auto param = params_free[idx_param];
            const auto found = _offsets_params.find(param);
            if (found == _offsets_params.end()) {
                throw std::runtime_error(param.get().str() + " not found in _offsets_params\n"
                                         + ERRMSG_PARAMS);
            }
            param_idx[idx_param] = (*found).second;
            if (not options) {
                const auto name = param.get().get_name();
                n_integral += name == IntegralParameterD::_name;
                n_frac += name == ProperFractionParameterD::_name;
            }
        }
 
        if (options) {
            const HessianOptions& opt = *options;
 
            auto loglike_grad = this->compute_loglike_grad(include_prior, false, false);
 
            size_t idx_param = 0;
            for (auto& paramref : params_free) {
                auto& param = paramref.get();
                const double value_init = param.get_value();
                const double value_transformed = param.get_value_transformed();
                // Make the finite differencing go away from the possible peak likelihood
                // This is probably a good idea even if return_negative isn't necessary
                double diff = std::copysign(std::abs(value_transformed * opt.findiff_frac) + opt.findiff_add,
                                            -loglike_grad[idx_param]);
                diff = finite_difference_param(param, diff);
                double diffabs = std::abs(diff);
                this->evaluate();
 
                auto loglike_grad2 = this->compute_loglike_grad(include_prior, false, false);
                for (size_t idx_col = 0; idx_col < n_params_free; ++idx_col) {
                    double dll_dx = loglike_grad2[idx_col] - loglike_grad[idx_col];
                    dll_dx *= opt.return_negative ? (1 - 2 * (dll_dx > 0)) / diffabs : 1 / diff;
                    hessian->set_value_unchecked(idx_param, idx_col, dll_dx);
                }
                param.set_value(value_init);
                idx_param++;
            }
        } else {
            if ((n_integral > 0) && (n_frac > 0)) {
                throw std::runtime_error(
                        "Cannot compute_hessian without options for model with free Integral and "
                        "ProperFraction "
                        "parameters (this is a bug); please supply options to use finite differencing.");
            }
 
            this->setup_evaluators(EvaluatorMode::jacobian);
            this->evaluate();
 
            hessian->fill(0);
 
            for (size_t idx_obs = 0; idx_obs < n_obs; ++idx_obs) {
                const auto& jacs = this->_grads.at(idx_obs);
                for (size_t idx_param1 = 0; idx_param1 < n_params_free; ++idx_param1) {
                    const auto& jac1 = jacs->at(param_idx[idx_param1]);
                    const size_t n_rows = jac1.get_n_rows();
                    const size_t n_cols = jac1.get_n_cols();
 
                    for (size_t idx_param2 = 0; idx_param2 < n_params_free; ++idx_param2) {
                        const auto& jac2 = jacs->at(param_idx[idx_param2]);
                        const size_t n_rows2 = jac2.get_n_rows();
                        const size_t n_cols2 = jac2.get_n_cols();
 
                        if ((n_rows != n_rows2) || (n_cols != n_cols2)) {
                            throw std::logic_error("n_rows,n_cols=" + to_string_float(n_rows) + ","
                                                   + to_string_float(n_cols)
                                                   + "!= n_rows2,n_cols2=" + to_string_float(n_rows2) + ","
                                                   + to_string_float(n_cols2));
                        }
                        double dx = 0;
                        for (size_t row = 0; row < n_rows; ++row) {
                            for (size_t col = 0; col < n_cols; ++col) {
                                // Note the sign convention of negative diagonal values
                                // at maximum likelihood
                                dx -= jac1.get_value(row, col) * jac2.get_value(row, col);
                            }
                        }
                        if (print) {
                            std::cout << "calling hessian->add_value(" << idx_param1 << "," << idx_param2
                                      << "," << dx << ")" << std::endl;
                        }
                        hessian->add_value(idx_param1, idx_param2, dx);
                    }
                }
            }
            if (include_prior && (_priors.size() > 0)) {
                std::map<ParamBaseCRef, size_t> param_idx_map = {};
                for (size_t idx_param = 0; idx_param < n_params_free; idx_param++) {
                    const auto param = params_free[idx_param];
                    if (_offsets_params.find(param) == _offsets_params.end()) {
                        throw std::runtime_error(param.get().str() + " not found in _offsets_params\n"
                                                 + ERRMSG_PARAMS);
                    }
                    param_idx_map[param] = idx_param;
                }
 
                size_t prior_size = 0;
                for (const auto& prior : _priors) {
                    prior_size += prior->size();
                }
                auto dll_prior = std::make_unique<Image>(n_params_free, prior_size);
                dll_prior->fill(0);
 
                size_t idx_prior = 0;
                for (const auto& prior : _priors) {
                    auto result = prior->evaluate(true);
 
                    for (const auto& [param_cref, values_base] : result.jacobians) {
                        size_t idx_residual = 0;
                        size_t idx_param = param_idx_map[param_cref];
                        for (const auto& value_base : values_base) {
                            dll_prior->add_value(idx_param, idx_prior + idx_residual++, value_base);
                        }
                    }
                    idx_prior += prior->size();
                }
 
                for (size_t row = 0; row < n_params_free; ++row) {
                    for (size_t row2 = 0; row2 < n_params_free; ++row2) {
                        double dll_dx = 0;
                        for (size_t col = 0; col < prior_size; ++col) {
                            dll_dx -= dll_prior->get_value(row, col) * dll_prior->get_value(row2, col);
                        }
                        hessian->add_value(row, row2, dll_dx);
                        if (print) {
                            std::cout << "calling hessian->add_value(" << row << "," << row2 << "," << dll_dx
                                      << ") for prior" << std::endl;
                        }
                    }
                }
            }
        }
 
        if (!transformed) {
            size_t idx_param = 0;
            for (auto& paramref : params_free) {
                auto& param = paramref.get();
                double value = param.get_value();
                param.set_transform(transforms[idx_param]);
                if (param.get_value() != value) {
                    throw std::logic_error("Param " + param.str() + " changed value from "
                                           + std::to_string(value) + " after dropping transform");
                }
                idx_param++;
            }
        }
 
        return hessian;
    }

compute_loglike_grad()

template<typename T, typename Image, typename Indices, typename Mask>

std::vector< double > lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::compute_loglike_grad ( bool include_prior = true, bool print = false, bool verify = false, double findiff_frac = 1e-5, double findiff_add = 1e-5, double rtol = 1e-3, double atol = 1e-8 )

inline

Compute the gradient (partial first derivative) of the log-likelihood for each free parameter.

Parameters

include_prior	Whether to include the prior likelihood(s) in the gradients.
print	Whether to print diagnostic/debugging information.
verify	Whether to verify the values by comparing to finite differences.
findiff_frac	The value of the finite difference increment, as a fraction of the parameter value.
findiff_add	The minimum value of the finite difference increment.
rtol	The allowed relative tolerance in the Jacobian as compared to the finite difference.
atol	The allowed absolute tolerance in the Jacobian as compared to the finite difference.

Returns

The gradient of the log likelihoods for the free parameters of the model, in the order returned by get_parameters.

Definition at line 704 of file model.h.

                                                                 {
        this->setup_evaluators(EvaluatorMode::loglike_grad, {}, {}, {}, nullptr, true, print);
        this->evaluate(print);
 
        auto filter_free = g2f::ParamFilter{false, true, true, true};
        auto params_free = this->get_parameters_const_new(&filter_free);
        params_free = nonconsecutive_unique<ParamBaseCRef>(params_free);
 
        const size_t n_params = this->_offsets_params.size();
        if (n_params != params_free.size()) {
            throw std::runtime_error("_offsets_params=" + str_map_ref<true>(_offsets_params)
                                     + ".size() = " + std::to_string(n_params)
                                     + " != params_free=" + str_iter_ref<true>(params_free) + ".size() + 1 ="
                                     + std::to_string(params_free.size()) + "; " + ERRMSG_PARAMS);
        }
        std::map<ParamBaseCRef, size_t> param_idx = {};
        for (size_t idx_param = 0; idx_param < n_params; idx_param++) {
            const auto param = params_free[idx_param];
            if (_offsets_params.find(param) == _offsets_params.end()) {
                throw std::runtime_error(param.get().str() + " not found in _offsets_params\n"
                                         + ERRMSG_PARAMS);
            }
            param_idx[param] = idx_param;
        }
        if (print) {
            std::cout << "params_free=" << str_iter_ref<true>(params_free) << std::endl;
            std::cout << "_offsets_params=" << str_map_ref<true>(_offsets_params) << std::endl;
            std::cout << "param_idx=" << str_map_ref<true>(param_idx) << std::endl;
        }
        std::vector<double> loglike_grads(n_params, 0.);
 
        size_t idx_obs = 0;
        for (const auto& grads_obs : _grads) {
            const auto& values = grads_obs->at(0);
            const size_t n_cols = values.get_n_cols();
            if (n_cols != n_params + 1) {
                throw std::logic_error(this->_data->at(idx_obs).get().str()
                                       + " loglike_grads n_cols=" + std::to_string(n_cols)
                                       + " != n_params (from offsets) = " + std::to_string(n_params + 1));
            }
            for (const auto& [paramref, col] : _offsets_params) {
                loglike_grads.at(param_idx[paramref]) += values.get_value(0, col);
            }
            idx_obs++;
        }
        if (include_prior) {
            for (const auto& prior : _priors) {
                auto result = prior->evaluate(true);
 
                // TODO: Confirm that this works with ShapePrior
                for (const auto& paramref : params_free) {
                    double dll_dx = result.compute_dloglike_dx(paramref);
                    loglike_grads.at(param_idx[paramref]) += dll_dx;
                }
            }
        }
        if (verify) {
            std::string errmsg;
            this->setup_evaluators(g2f::EvaluatorMode::loglike);
            auto loglike = this->evaluate();
 
            auto params = this->get_parameters_new(&filter_free);
            params = nonconsecutive_unique<ParamBaseRef>(params);
 
            for (size_t idx_param = 0; idx_param < n_params; idx_param++) {
                auto& param = params[idx_param].get();
                double value = param.get_value_transformed();
                double diff = std::copysign(std::abs(value * findiff_frac) + findiff_add,
                                            loglike_grads[idx_param]);
                diff = finite_difference_param(param, diff);
                std::vector<double> loglike_new_plus;
                std::vector<double> loglike_new_minus;
                try {
                    param.set_value_transformed(value + diff / 2.);
                    loglike_new_plus = this->evaluate();
                    param.set_value_transformed(value - diff / 2.);
                    loglike_new_minus = this->evaluate();
                } catch (std::runtime_error& e) {
                    param.set_value_transformed(value + diff);
                    loglike_new_plus = this->evaluate();
                    loglike_new_minus = loglike;
                }
                auto dloglike = (sum_iter(loglike_new_plus) - sum_iter(loglike_new_minus)) / diff;
                auto close = isclose(dloglike, loglike_grads[idx_param], rtol, atol);
                param.set_value_transformed(value);
                if (!close.isclose) {
                    double dll_dx_exact = loglike_grads[idx_param];
                    double dlp_dx_sum = 0;
                    for (const auto& prior : _priors) {
                        double dlp_dx = prior->evaluate(true).compute_dloglike_dx(param, true);
                        dlp_dx_sum += dlp_dx;
                    }
                    double dlp_dx_findiff = (loglike_new_plus.back() - loglike_new_minus.back()) / diff;
                    errmsg += param.str() + " failed loglike_grad verification; isclose=" + close.str()
                              + " from findiff=" + to_string_float(dloglike) + " vs "
                              + to_string_float(dll_dx_exact)
                              + " (diff = " + to_string_float(dll_dx_exact - dloglike)
                              + " , ratio = " + to_string_float(dll_dx_exact / dloglike)
                              + "); dll_dx_prior=" + to_string_float(dlp_dx_sum)
                              + " vs findiff: " + to_string_float(dlp_dx_findiff) + "\n";
                }
            }
            if (errmsg.size() > 0) throw std::logic_error(errmsg);
        }
 
        return loglike_grads;
    }

evaluate()

template<typename T, typename Image, typename Indices, typename Mask>

std::vector< double > lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::evaluate ( bool print = false, bool normalize_loglike = false )

inline

Evaluate the model for every Observation in _data.

Parameters

print	Whether to print diagnostic statements to stdout.
normalize_loglike	Whether to include the normalizing (variance-dependent) term in the log likelihood. If false, the log likelihood for a model with no residuals is 0.

Returns

The log likelihood of each observation, followed by the summed log likelihood of the priors.

Definition at line 1019 of file model.h.

                                                                                   {
        if (!_is_setup) throw std::runtime_error("Can't call evaluate before setup_evaluators");
 
        std::vector<double> result(_size + 1);
        bool is_loglike_grad = this->_mode == EvaluatorMode::loglike_grad;
 
        if (is_loglike_grad) {
            // loglike_grad adds values, and they have to be reset to zero first
            // This must be done for all observations in this case, because
            // most parameters change grads in all observations
            for (size_t idx = 0; idx < _size; ++idx) {
                this->_grads[idx]->at(0).fill(0);
            }
        }
 
        for (size_t idx = 0; idx < _size; ++idx) {
            result[idx] = this->_evaluate_observation(idx, print, is_loglike_grad);
        }
        if ((this->_mode == EvaluatorMode::loglike) || is_loglike_grad
            || (this->_mode == EvaluatorMode::loglike_image) || (this->_mode == EvaluatorMode::jacobian)) {
            result[_size] = this->_evaluate_priors(print, normalize_loglike);
        }
 
        return result;
    }

evaluate_observation()

template<typename T, typename Image, typename Indices, typename Mask>

double lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::evaluate_observation ( size_t idx )

inline

Evaluate a single observation with the given index in _data.

Parameters

idx	The numeric index of the observation.

Returns

The log likelihood of the model for that observation.

Definition at line 1051 of file model.h.

                                            {
        _check_obs_idx(idx);
        return _evaluate_observation(_evaluators[idx]);
    }

get_data()

template<typename T, typename Image, typename Indices, typename Mask>

std::shared_ptr< const ModelData > lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_data ( ) const

inline

Return _data.

Definition at line 1057 of file model.h.

{ return _data; }

get_gaussians()

template<typename T, typename Image, typename Indices, typename Mask>

std::unique_ptr< const lsst::gauss2d::Gaussians > lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_gaussians ( const Channel & channel ) const

inlineoverridevirtual

Return the vector of Gaussian sub-components controlled by this model.

Parameters

channel

The Channel to return Gaussians for.

Returns

The Gaussians controlled by this model.

Implements lsst::gauss2d::fit::ParametricModel.

Definition at line 1061 of file model.h.

                                                                                                   {
        std::vector<std::optional<const lsst::gauss2d::Gaussians::Data>> in;
        // TODO: Rethink this; it's not sufficient unless sources are single component gaussians
        in.reserve(_sources.size());
 
        for (auto& source : _sources) {
            auto data = source->get_gaussians(channel)->get_data();
            in.emplace_back(data);
        }
 
        return std::make_unique<lsst::gauss2d::Gaussians>(in);
    }

get_loglike_const_terms()

template<typename T, typename Image, typename Indices, typename Mask>

std::vector< double > lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_loglike_const_terms ( )

inline

Get the constant (variance-dependent) terms of the log likelihood for each observation.

Returns

A vector with the constant terms for each observation, and a final term for the priors.

Definition at line 1079 of file model.h.

                                                {
        const size_t n_data = this->size();
        if (this->_likelihood_const_terms.empty()) {
            this->_likelihood_const_terms.resize(n_data + 1);
 
            for (size_t idx = 0; idx < n_data; ++idx) {
                const auto& sigma_inv = this->_data[idx].get_sigma_inverse();
                double loglike = 0;
                const size_t n_rows = sigma_inv.get_n_rows();
                const size_t n_cols = sigma_inv.get_n_cols();
                for (size_t col = 0; col < n_cols; ++col) {
                    for (size_t row = 0; row < n_rows; ++row) {
                        loglike += LOG_1 - log(sigma_inv._get_value(row, col));
                    }
                }
                _likelihood_const_terms[idx] = loglike;
            }
        }
        double loglike = 0;
        for (const auto& prior : this->_priors) {
            auto values = prior->get_loglike_const_terms();
            for (const auto value : values) loglike += value;
        }
        _likelihood_const_terms[n_data] = loglike;
        return _likelihood_const_terms;
    }

get_mode()

template<typename T, typename Image, typename Indices, typename Mask>

EvaluatorMode lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_mode ( ) const

inline

Definition at line 1059 of file model.h.

{ return _mode; }

get_n_gaussians()

template<typename T, typename Image, typename Indices, typename Mask>

size_t lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_n_gaussians ( const Channel & channel ) const

inlineoverridevirtual

Return the number of Gaussian sub-components controlled by this model.

Implements lsst::gauss2d::fit::ParametricModel.

Definition at line 1106 of file model.h.

                                                                  {
        size_t n_g = 0;
        for (auto& source : _sources) n_g += source->get_n_gaussians(channel);
        return n_g;
    }

get_offsets_parameters()

template<typename T, typename Image, typename Indices, typename Mask>

std::vector< std::pair< ParamBaseCRef, size_t > > lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_offsets_parameters ( ) const

inline

Definition at line 1115 of file model.h.

                                                                             {
        std::vector<std::pair<ParamBaseCRef, size_t>> offsets = {};
        for (auto& [param, offset] : _offsets_params) {
            offsets.emplace_back(param, offset);
        }
        return offsets;
    }

get_outputs()

template<typename T, typename Image, typename Indices, typename Mask>

std::vector< std::shared_ptr< Image > > lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_outputs ( ) const

inline

Return _outputs (output Image instances for each Observation in _data)

Definition at line 1113 of file model.h.

{ return _outputs; }

get_parameters()

template<typename T, typename Image, typename Indices, typename Mask>

ParamRefs & lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_parameters ( ParamRefs & params, ParamFilter * filter = nullptr ) const

inlineoverridevirtual

Add Parameter refs matching the filter to a vector, in order.

Parameters

params	The vector to add to.
filter	The filter to apply to this Object's parameters.

Returns

A ref to params (for method chaining)

Implements lsst::gauss2d::fit::Parametric.

Definition at line 1123 of file model.h.

                                                                                               {
        for (auto& source : _sources) source->get_parameters(params, filter);
        for (auto& psfmodel : _psfmodels) psfmodel->get_parameters(params, filter);
        _data->get_parameters(params, filter);
        return params;
    }

get_parameters_const()

template<typename T, typename Image, typename Indices, typename Mask>

ParamCRefs & lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_parameters_const ( ParamCRefs & params, ParamFilter * filter = nullptr ) const

inlineoverridevirtual

Same as get_parameters(), but for const refs.

Implements lsst::gauss2d::fit::Parametric.

Definition at line 1130 of file model.h.

                                                                                                       {
        for (auto& source : _sources) source->get_parameters_const(params, filter);
        for (auto& psfmodel : _psfmodels) psfmodel->get_parameters_const(params, filter);
        _data->get_parameters_const(params, filter);
        return params;
    }

get_parameters_const_new()

ParamCRefs lsst::gauss2d::fit::Parametric::get_parameters_const_new ( ParamFilter * filter = nullptr ) const

inlineinherited

Same as get_parameters_const(), but returning a new vector.

Definition at line 33 of file parametric.h.

                                                                             {
        ParamCRefs params{};
        get_parameters_const(params, filter);
        return params;
    }

get_parameters_new()

ParamRefs lsst::gauss2d::fit::Parametric::get_parameters_new ( ParamFilter * filter = nullptr ) const

inlineinherited

Same as get_parameters(), but returning a new vector.

Definition at line 27 of file parametric.h.

                                                                      {
        ParamRefs params{};
        get_parameters(params, filter);
        return params;
    }

get_parameters_observation()

template<typename T, typename Image, typename Indices, typename Mask>

ParamRefs & lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_parameters_observation ( ParamRefs & params, size_t idx, ParamFilter * filter = nullptr ) const

inline

Same as get_parameters(), but for a single Observation with index idx in _data.

Definition at line 1138 of file model.h.

                                                                               {
        _check_obs_idx(idx);
        for (auto& source : _sources) source->get_parameters(params, filter);
        _psfmodels[idx]->get_parameters(params, filter);
        _data->at(idx).get().get_parameters(params, filter);
        return params;
    }

get_parameters_observation_const()

template<typename T, typename Image, typename Indices, typename Mask>

ParamCRefs & lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_parameters_observation_const ( ParamCRefs & params, size_t idx, ParamFilter * filter = nullptr ) const

inline

Same as get_parameters_const(), but for a single Observation with index idx in _data.

Definition at line 1148 of file model.h.

                                                                                      {
        _check_obs_idx(idx);
        for (auto& source : _sources) source->get_parameters_const(params, filter);
        _psfmodels[idx]->get_parameters_const(params, filter);
        _data->at(idx).get().get_parameters_const(params, filter);
        return params;
    }

get_priors()

template<typename T, typename Image, typename Indices, typename Mask>

Priors lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_priors ( ) const

inline

Return _priors, the list of Prior instances.

Definition at line 1158 of file model.h.

{ return _priors; }

get_psfmodels()

template<typename T, typename Image, typename Indices, typename Mask>

PsfModels lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_psfmodels ( ) const

inline

Return _psfmodels, the list of PsfModel instances for each Observation in _data.

Definition at line 1161 of file model.h.

{ return _psfmodels; }

get_sources()

template<typename T, typename Image, typename Indices, typename Mask>

Sources lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::get_sources ( ) const

inline

Return _sources, the list of Source instances for each Observation in _data.

Definition at line 1164 of file model.h.

{ return _sources; }

null_str()

static std::string_view lsst::gauss2d::Object::null_str ( const std::string_view & namespace_separator )

inlinestaticinherited

Definition at line 49 of file object.h.

                                                                              {
        return namespace_separator == CC_NAMESPACE_SEPARATOR ? "nullptr" : NULL_STR_GENERAL;
    }

repr()

template<typename T, typename Image, typename Indices, typename Mask>

std::string lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::repr ( bool name_keywords = false, std::string_view namespace_separator = Object::CC_NAMESPACE_SEPARATOR ) const

inlineoverridevirtual

Return a full, callable string representation of this.

Parameters

name_keywords	Whether to prefix arguments with "{name}=", where name is the arg name in the header (as with keyword arguments in Python).
namespace_separator	The string to use to delimit namespaces, i.e. :: in C++ and . in Python.

Returns

A callable string representation of this, which should return an an identical object to this.

Note

The representation with name_keywords=false must be callable in C++. The representation with name_keywords=true should be callable in Python, if there are any bindings.

Implements lsst::gauss2d::Object.

Definition at line 1377 of file model.h.

                                                                                                       {
        std::string str = type_name_str<Model>(false, namespace_separator) + "("
                          + (name_keywords ? "data=" : "") + _data->repr(name_keywords, namespace_separator)
                          + ", " + (name_keywords ? "psfmodels=[" : "[");
        for (const auto& x : _psfmodels) str += x->repr(name_keywords, namespace_separator) + ",";
        str += (name_keywords ? "], sources=[" : "], [");
        for (const auto& s : _sources) str += s->repr(name_keywords, namespace_separator) + ",";
        str += (name_keywords ? "], priors=[" : "], [");
        return str + "])";
    }

set_extra_param_factors()

template<typename T, typename Image, typename Indices, typename Mask>

void lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::set_extra_param_factors ( const Channel & channel, ExtraParamFactors & factors, size_t index ) const

inlineoverridevirtual

Set extra Parameter gradient factors in an existing map.

Parameters

channel	The Channel to set factors for.
factors	The ExtraParamFactors to set factors for.
index	The index to begin setting factors at.

Implements lsst::gauss2d::fit::ParametricModel.

Definition at line 1166 of file model.h.

                                                              {
        for (auto& source : _sources) {
            source->set_extra_param_factors(channel, factors, index);
            index += source->get_n_gaussians(channel);
        }
    }

set_grad_param_factors()

template<typename T, typename Image, typename Indices, typename Mask>

void lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::set_grad_param_factors ( const Channel & channel, GradParamFactors & factors, size_t index ) const

inlineoverridevirtual

Set Parameter gradient factors in an existing map.

Parameters

channel	The Channel to set factors for.
factors	The GradParamFactors to set factors for.
index	The index to begin setting factors at.

Implements lsst::gauss2d::fit::ParametricModel.

Definition at line 1174 of file model.h.

                                                             {
        for (auto& source : _sources) {
            source->set_grad_param_factors(channel, factors, index);
            index += source->get_n_gaussians(channel);
        }
    }

setup_evaluators()

template<typename T, typename Image, typename Indices, typename Mask>

void lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::setup_evaluators ( EvaluatorMode mode = EvaluatorMode::image, std::vector< std::vector< std::shared_ptr< Image > > > outputs = {}, std::vector< std::shared_ptr< Image > > residuals = {}, std::vector< std::shared_ptr< Image > > outputs_prior = {}, std::shared_ptr< Image > residuals_prior = nullptr, bool force = false, bool print = false )

inline

Setup Evaluator instances for every Observation in _data using the given EvaluatorMode.

Parameters

mode	The EvaluatorMode to use for all Evaluator instances.
outputs	A vector of vectors of Image outputs for each Evaluator (created if empty and needed).
residuals	A vector of residual Images for each Evaluator (created if empty and needed).
outputs_prior	A vector of prior output (Jacobian) Images for each Evaluator (created if empty and needed).
residuals_prior	A vector of prior residual Images for each Evaluator (created if empty and needed).
force	Whether to force setting up even if already set up in the same mode
print	Whether to print diagnostic statements to stdout.

Note

Different modes require different sized vectors for outputs EvaluatorMode::jacobian requires one Image per free parameter per Observation. EvaluatorMode::loglike_grad requires only one Image(n_rows=1, n_cols=n_params_free)

Definition at line 1199 of file model.h.

                                                                               {},
                          std::vector<std::shared_ptr<Image>> residuals = {},
                          std::vector<std::shared_ptr<Image>> outputs_prior = {},
                          std::shared_ptr<Image> residuals_prior = nullptr, bool force = false,
                          bool print = false) {
        const size_t n_outputs = outputs.size();
        const bool has_outputs = n_outputs > 0;
        if (has_outputs) {
            if (n_outputs != size()) {
                throw std::invalid_argument("outputs.size()=" + std::to_string(n_outputs)
                                            + "!=this->size()=" + std::to_string(size()));
            }
        }
        const size_t n_residuals = residuals.size();
        const bool has_residuals = n_residuals > 0;
        if (has_residuals) {
            if (n_residuals != size()) {
                throw std::invalid_argument("residuals.size()=" + std::to_string(n_residuals)
                                            + "!=this->size()=" + std::to_string(size()));
            }
        }
        if (force || (!_is_setup || (mode != _mode))) {
            _evaluators.clear();
            _grads.clear();
            _outputs.clear();
            _outputs_prior.clear();
            _residuals_prior = nullptr;
            _evaluators.reserve(size());
            _is_setup = false;
 
            const bool is_jacobian = mode == EvaluatorMode::jacobian;
            const bool is_loglike_grad = mode == EvaluatorMode::loglike_grad;
 
            std::vector<ParamBaseCRef> params_free = {};
            if (is_jacobian || is_loglike_grad) {
                auto filter_free = g2f::ParamFilter{false, true, true, true};
                this->get_parameters_const(params_free, &filter_free);
                params_free = nonconsecutive_unique<ParamBaseCRef>(params_free);
                _n_params_free = params_free.size();
                if (print) {
                    std::cout << "n_params_free=" << _n_params_free << " from params_free:" << std::endl;
                    std::cout << str_iter_ref<true>(params_free) << std::endl;
                }
 
                _grads.reserve(size());
                _offsets_params.clear();
            }
 
            size_t n_residuals_prior = 0;
            bool do_residuals_prior = (_size_priors > 0)
                                      && (is_loglike_grad || is_jacobian || (mode == EvaluatorMode::loglike)
                                          || (_mode == EvaluatorMode::loglike_image));
            if (do_residuals_prior) {
                for (const auto& prior : this->_priors) {
                    n_residuals_prior += prior->size();
                }
 
                _residuals_prior
                        = residuals_prior == nullptr
                                  ? (_residuals_prior = std::make_shared<Image>(1, n_residuals_prior))
                                  : std::move(residuals_prior);
                if ((_residuals_prior->get_n_rows() != 1)
                    || (_residuals_prior->get_n_cols() != n_residuals_prior)) {
                    throw std::invalid_argument("residuals_prior=" + _residuals_prior->str()
                                                + " rows,cols != 1," + std::to_string(n_residuals_prior));
                }
                if (print) std::cout << "residuals_prior built/validated" << std::endl;
                for (size_t col = 0; col < n_residuals_prior; ++col) {
                    _residuals_prior->set_value(0, col, 0);
                }
                if (print) std::cout << "residuals_prior reset" << std::endl;
 
                if (is_jacobian) {
                    if (print) {
                        std::cout << "setup_evaluators jacobian called for model:" << std::endl;
                        std::cout << this->str() << std::endl;
                    }
 
                    const size_t n_params_jac = _n_params_free + 1;
                    if (outputs_prior.size() == 0) {
                        for (size_t idx_jac = 0; idx_jac < n_params_jac; ++idx_jac) {
                            _outputs_prior.emplace_back(std::make_shared<Image>(1, n_residuals_prior));
                        }
                        if (print) std::cout << "outputs_prior built" << std::endl;
                    } else if (outputs_prior.size() != n_params_jac) {
                        throw std::invalid_argument(
                                "jacobian outputs_prior->size()=" + std::to_string(outputs_prior.size())
                                + " != (n_params_free + 1 = " + std::to_string(n_params_jac) + ")");
                    }
                    _outputs_prior.reserve(n_params_jac);
                    size_t idx_jac = 0;
                    for (auto& output_prior : outputs_prior) {
                        if (output_prior == nullptr) {
                            throw std::invalid_argument("output_prior[" + std::to_string(idx_jac)
                                                        + "] is null");
                        }
                        if ((output_prior->get_n_rows() != 1)
                            || (output_prior->get_n_cols() != n_residuals_prior)) {
                            throw std::invalid_argument("outputs_prior[" + std::to_string(idx_jac)
                                                        + "]=" + output_prior->str() + " rows,cols != 1,"
                                                        + std::to_string(n_residuals_prior));
                        }
                        for (size_t col = 0; col < n_residuals_prior; ++col) {
                            output_prior->set_value(0, col, 0);
                        }
                        _outputs_prior.emplace_back(std::move(output_prior));
                        idx_jac++;
                    }
                    if (print) std::cout << "outputs_prior validated" << std::endl;
                }
            }
            if (!is_jacobian) {
                _outputs.reserve(size());
            }
 
            const auto& channels = _data->get_channels();
 
            if (print) std::cout << "making evaluators" << std::endl;
 
            for (size_t idx = 0; idx < _size; ++idx) {
                auto result = this->_make_evaluator(idx, mode, outputs,
                                                    has_residuals ? residuals[idx] : nullptr, print);
                _evaluators.emplace_back(std::move(result.first));
                if (is_jacobian || is_loglike_grad) {
                    _grads.emplace_back(std::move(result.second.second));
                } else {
                    _outputs.emplace_back(std::move(result.second.first));
                }
            }
 
            if (print) std::cout << "evaluators made" << std::endl;
 
            std::set<ParamBaseCRef> params_prior;
            size_t idx_prior = 0;
 
            // Check that the priors are reasonable.
            for (const auto& prior : _priors) {
                auto eval_prior = prior->evaluate(is_jacobian);
                size_t n_resid = prior->size();
                size_t n_resid_eval = eval_prior.residuals.size();
                if (n_resid != eval_prior.residuals.size()) {
                    throw std::logic_error(prior->str() + ".size()=" + std::to_string(n_resid)
                                           + " != evaluate(true).residuals.size()="
                                           + std::to_string(n_resid_eval));
                }
                if (is_jacobian) {
                    size_t idx_param = 0;
                    for (const auto& param_jacob : eval_prior.jacobians) {
                        const auto& param = param_jacob.first;
                        if (params_prior.find(param) != params_prior.end()) {
                            throw std::invalid_argument(
                                    "_priors[" + std::to_string(idx_prior) + "]=" + prior->str()
                                    + ".evaluate(true)[" + std::to_string(idx_param)
                                    + "]=" + param.get().str() + " already in prior param set."
                                    + " Did you apply multiple priors to the same parameter?");
                        }
                        idx_param++;
                    }
                }
                idx_prior++;
            }
            if (print && (_priors.size() > 0)) std::cout << "priors validated" << std::endl;
            if (_size_priors > 0) {
                if ((_residuals_prior != nullptr) && (_residuals_prior->get_n_cols() != n_residuals_prior)) {
                    throw std::invalid_argument("residuals_prior=" + _residuals_prior->str() + " n_residuals="
                                                + std::to_string(n_residuals_prior) + "!= get_n_cols()="
                                                + std::to_string(_residuals_prior->get_n_cols()));
                }
            }
            _is_setup = true;
            _mode = mode;
        }
    }

size()

template<typename T, typename Image, typename Indices, typename Mask>

size_t lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::size ( ) const

inline

Get the size of this->_data.

Definition at line 1375 of file model.h.

{ return _size; }

str()

template<typename T, typename Image, typename Indices, typename Mask>

std::string lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::str ( ) const

inlineoverridevirtual

Return a brief, human-readable string representation of this.

Implements lsst::gauss2d::Object.

Definition at line 1389 of file model.h.

                                   {
        std::string str = type_name_str<Model>(true) + "(data=" + _data->str() + ", psfmodels=[";
        for (const auto& x : _psfmodels) str += x->str() + ",";
        str += "], sources=[";
        for (const auto& x : _sources) str += x->str() + ",";
        str += "], priors=[";
        for (const auto& x : _priors) str += x->str() + ",";
        return str + "])";
    }

verify_jacobian()

template<typename T, typename Image, typename Indices, typename Mask>

std::vector< std::string > lsst::gauss2d::fit::Model< T, Image, Indices, Mask >::verify_jacobian ( double findiff_frac = 1e-5, double findiff_add = 1e-5, double rtol = 1e-3, double atol = 1e-3, double max_ll_diff = 0 )

inline

Verify that the Jacobian is correct by comparing to finite differences.

Parameters

findiff_frac	The value of the finite difference increment, as a fraction of the parameter value.
findiff_add	The minimum value of the finite difference increment.
rtol	The allowed relative tolerance in the Jacobian as compared to the finite difference.
atol	The allowed absolute tolerance in the Jacobian as compared to the finite difference.
max_ll_diff	The maximum allowed difference between equivalent log-likelihood evaluations. Must be >= 0 and should only be a few orders of magnitude larger than machine epsilon.

Returns

A vector of error messages, one for each Parameter that failed verification.

Note

Verification is done using isclose(), which is modelled after Python's numpy.isclose.

Definition at line 1412 of file model.h.

                                                                                                             {
        if (!(max_ll_diff >= 0)) {
            throw std::invalid_argument("max_ll_diff=" + to_string_float(max_ll_diff) + " !> 0");
        }
        if (_mode != EvaluatorMode::jacobian) {
            this->setup_evaluators(EvaluatorMode::jacobian);
        }
 
        const size_t n_obs = this->size();
        ParamFilter filter{false, true};
        std::vector<std::string> errors;
 
        for (size_t idx = 0; idx < n_obs; ++idx) {
            const auto& grads = *(_grads.at(idx));
            const auto& observation = _data->at(idx).get();
            const auto& sigma_inv = observation.get_sigma_inverse();
 
            const Channel& channel = observation.get_channel();
            const size_t n_cols = observation.get_n_cols();
            const size_t n_rows = observation.get_n_rows();
 
            filter.channel = channel;
 
            ParamRefs params{};
            this->get_parameters_observation(params, idx, &filter);
            params = nonconsecutive_unique(params);
 
            size_t idx_param_max = 0;
            for (const auto& param : params) {
                auto found = _offsets_params.find(param);
                if (found == _offsets_params.end()) {
                    throw std::runtime_error(
                            param.get().str()
                            + " not found in offsets; was it freed after setup_evaluators was called?");
                }
                size_t idx_param = found->second;
                if (idx_param > idx_param_max) idx_param_max = idx_param;
            }
 
            if (!(idx_param_max < grads.size())) {
                throw std::runtime_error("idx_param_max=" + std::to_string(idx_param_max)
                                         + " !< grads.size()=" + std::to_string(grads.size())
                                         + "; is the jacobian array large enough?");
            }
 
            double loglike_jac = this->_evaluate_observation(idx);
            const auto& grad = grads[0];
            size_t n_nonzero = 0;
            for (unsigned int i = 0; i < n_cols; ++i) {
                for (unsigned int j = 0; j < n_rows; ++j) {
                    n_nonzero += grad.get_value(j, i) != 0;
                }
            }
            if (n_nonzero > 0) {
                errors.push_back("n_nonzero grads[0] entries=" + std::to_string(n_nonzero));
            }
 
            auto result = _make_evaluator(idx, EvaluatorMode::loglike_image);
            double loglike_img = result.first->loglike_pixel();
            auto is_close_ll = isclose(loglike_img, loglike_jac, 0., max_ll_diff);
            if (!is_close_ll.isclose) {
                errors.push_back("loglike_jac=" + to_string_float(loglike_jac) + ", loglike_img="
                                 + to_string_float(loglike_img) + " !close (isclose=" + is_close_ll.str()
                                 + ") for obs[" + std::to_string(idx) + "]:");
            }
 
            const auto& image_current = *(result.second.first);
 
            for (auto& param_ref : params) {
                auto& param = param_ref.get();
                const size_t idx_param = _offsets_params.find(param_ref)->second;
                const auto& grad = grads[idx_param];
 
                const double value_init = param.get_value();
                const double value_transformed = param.get_value_transformed();
                double diff = value_transformed * findiff_frac;
                if (std::abs(diff) < findiff_add) diff = findiff_add;
                diff = finite_difference_param(param, diff);
 
                auto result_diff = _make_evaluator(idx, EvaluatorMode::image);
 
                result_diff.first->loglike_pixel();
 
                size_t n_failed = 0;
                std::vector<double> ratios = {};
                for (unsigned int i = 0; i < n_cols; ++i) {
                    for (unsigned int j = 0; j < n_rows; ++j) {
                        double delta = sigma_inv.get_value(j, i) / diff
                                       * (result_diff.second.first->get_value(j, i)
                                          - image_current.get_value(j, i));
                        double jac = grad.get_value(j, i);
 
                        auto close = isclose(jac, delta, rtol, atol);
                        if (!close.isclose) {
                            n_failed++;
                            ratios.push_back(jac / delta);
                        }
                    }
                }
 
                param.set_value(value_init);
                const double value_new = param.get_value();
                if (value_new != value_init) {
                    throw std::logic_error("Could not return param=" + param.str()
                                           + " to value_init=" + to_string_float(value_init) + "; diff="
                                           + to_string_float(value_new - value_init) + "); check limits");
                }
                if (n_failed > 0) {
                    std::sort(ratios.begin(), ratios.end());
                    double median = ratios[ratios.size() / 2];
                    errors.push_back("Param[" + std::to_string(idx_param) + "]=" + param.str()
                                     + " failed for obs[" + std::to_string(idx)
                                     + "]: n_failed=" + std::to_string(n_failed)
                                     + "; median evaluated/expected=" + std::to_string(median));
                }
            }
        }
        filter.channel = std::nullopt;
 
        ParamRefs params{};
        this->get_parameters_observation(params, 0, &filter);
        params = nonconsecutive_unique(params);
 
        for (const auto& prior : this->_priors) {
            auto result_base = prior->evaluate(true);
            auto residuals_base = result_base.residuals;
            const size_t n_values = residuals_base.size();
 
            for (const auto& [param_cref, values_base] : result_base.jacobians) {
                if (values_base.size() != n_values) {
                    throw std::logic_error("values_base.size()=" + std::to_string(values_base.size())
                                           + " != residuals_base.size()=" + std::to_string(n_values));
                }
 
                // TODO: Determine if it's possible to get mutable ref from set.find
                // (set.find only ever seems to return cref, so maybe not)
                auto found = std::find(params.begin(), params.end(), param_cref);
                if (found == params.end()) {
                    throw std::logic_error("Couldn't find " + param_cref.get().str() + " in all params");
                }
                auto& param = (*found).get();
 
                const double value_init = param.get_value();
                const double value_transformed = param.get_value_transformed();
                double delta = value_transformed * findiff_frac;
                if (std::abs(delta) < findiff_add) delta = findiff_add;
                delta = finite_difference_param(param, delta);
                auto result_new = prior->evaluate(true);
 
                const size_t idx_param = _offsets_params.at(param);
 
                std::vector<double> ratios = {};
                size_t n_failed = 0;
                size_t idx_value = 0;
                for (const double jac_base : values_base) {
                    double jac_findiff
                            = (result_new.residuals[idx_value] - result_base.residuals[idx_value]) / delta;
                    auto close = isclose(jac_base, jac_findiff, rtol, atol);
                    if (!close.isclose) {
                        n_failed++;
                        prior->evaluate(true);
                        ratios.push_back(jac_base / jac_findiff);
                    }
                    idx_value++;
                }
                param.set_value(value_init);
                const double value_new = param.get_value();
                if (value_new != value_init) {
                    throw std::logic_error("Could not return param=" + param.str()
                                           + " to value_init=" + to_string_float(value_init) + "; diff="
                                           + to_string_float(value_new - value_init) + "); check limits");
                }
                if (n_failed > 0) {
                    std::sort(ratios.begin(), ratios.end());
                    double median = ratios[ratios.size() / 2];
                    errors.push_back("Param[" + std::to_string(idx_param) + "]=" + param.str()
                                     + " failed for prior=" + prior->str()
                                     + ": n_failed=" + std::to_string(n_failed)
                                     + "; median evaluated/expected=" + std::to_string(median));
                }
            }
        }
        return errors;
    }

Member Data Documentation

CC_NAMESPACE_SEPARATOR

std::string_view lsst::gauss2d::Object::CC_NAMESPACE_SEPARATOR = "::"

staticconstexprinherited

The C++ namespace separator.

Definition at line 45 of file object.h.

NULL_STR_GENERAL

std::string_view lsst::gauss2d::Object::NULL_STR_GENERAL = "None"

staticconstexprinherited

Definition at line 46 of file object.h.

PY_NAMESPACE_SEPARATOR

std::string_view lsst::gauss2d::Object::PY_NAMESPACE_SEPARATOR = "."

staticconstexprinherited

Definition at line 47 of file object.h.

---

The documentation for this class was generated from the following file:

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-10.0.0/Linux64/gauss2d_fit/gb58c049af0+f03b321e39/include/lsst/gauss2d/fit/model.h