Evaluate the log likelihood and residual-dependent terms.
127 {
129 this->_ellipse->get_rho());
131
132 double size_maj = ellipse_major.get_r_major();
133 double axrat = size_maj == 0 ? 0 : ellipse_major.get_axrat();
134
135 if (!(axrat >= 0)) {
137 }
138
139 auto result = PriorEvaluation(0, {}, {},
false);
140
141 if (this->_prior_size != nullptr) {
142 double size_maj_floor = this->_options->get_size_maj_floor();
143 size_maj = sqrt(size_maj * size_maj * axrat + size_maj_floor * size_maj_floor);
144 size_maj = _prior_size->get_mean_parameter().get_transform().forward(size_maj);
145 double sigma = _prior_size->get_stddev();
146 double residual = (size_maj - _prior_size->get_mean_parameter().get_value_transformed()) /
sigma;
150 + ", _prior_size=" + _prior_size->str());
151 }
152 result.residuals.emplace_back(residual);
153 result.loglike += normalize ?
logpdf_norm(residual, 1.0) : -residual * residual / 2.;
154 }
155
156 if (this->_prior_axrat != nullptr) {
157 double axrat_floor = this->_options->get_axrat_floor();
158 axrat = sqrt(axrat_floor * axrat_floor + axrat * axrat);
159 if (axrat > 1) {
160 axrat = 1;
161 } else if (!(axrat >= 0)) {
164 }
165 const auto& transform_param = _prior_axrat->get_mean_parameter().get_transform();
166 axrat = transform_param.forward(axrat);
167 double sigma = _prior_axrat->get_stddev();
168 double residual = (axrat - _prior_axrat->get_mean_parameter().get_value_transformed()) /
sigma;
170 auto str = this->
str();
173 + ", _prior_axrat=" + _prior_axrat->str());
174 }
175
176
177 result.residuals.emplace_back(residual);
178 result.loglike += normalize ?
logpdf_norm(residual, 1.0) : -residual * residual / 2.;
179 }
180
181 if (calc_jacobians) {
183 for (auto& paramref : params) {
184 auto& param = paramref.get();
185 double value_init = param.get_value();
186 double value_trans = param.get_value_transformed();
188 double value_new = param.get_value_transformed();
189
192 try {
193 param.set_value_transformed(value_trans + delta / 2.);
195 param.set_value_transformed(value_trans - delta / 2.);
198 param.set_value_transformed(value_new);
200 residuals_old =
result.residuals;
201 }
202
203 param.set_value(value_init);
204 value_new = param.get_value();
205 if (value_new != value_init) {
209 }
211 for (
size_t idx = 0; idx <
result.residuals.size(); ++idx) {
212 jacobians[idx] = (residuals_new[idx] - residuals_old[idx]) / delta;
213 }
214 result.jacobians[paramref] = jacobians;
215 }
216 }
217
219}
afw::table::Key< double > sigma
An Ellipse with sigma_x, sigma_y, and rho values.
An Ellipse with r_major, axrat and angle values.
std::vector< double > residuals
std::vector< ParamBaseRef > ParamRefs
double finite_difference_param(ParamBase ¶m, double delta)
double logpdf_norm(T residual, T sigma)
std::vector< T > nonconsecutive_unique(const std::vector< T > &vec)
std::string to_string_float(const T value, const int precision=6, const bool scientific=true)