97 _size(_data == nullptr ? 0 : _data->
size()),
98 _size_priors(priors.
size()) {
106 _outputs.resize(
size);
109 for (
auto& psfmodel : psfmodels) {
110 if (psfmodel ==
nullptr)
112 _psfcomps.
push_back(psfmodel->get_gaussians());
117 _priors.reserve(_size_priors);
119 for (
auto& prior : priors) {
120 if (prior ==
nullptr)
122 _priors.push_back(prior);
126 _sources.
reserve(sources.size());
128 for (
auto& source : sources) {
129 if (source ==
nullptr)
135 _gaussians_srcs = this->_get_gaussians_srcs();
136 _factors_extra_in.resize(_size);
137 _factors_grad_in.resize(_size);
138 _map_extra_in.
resize(_size);
139 _map_grad_in.
resize(_size);
151 bool _is_setup =
false;
156 size_t _n_params_free = 0;
168 void _check_obs_idx(
size_t idx)
const {
169 if (!(idx < _size)) {
178 size_t n_gaussians_src = 0;
180 for (
const Channel& channel : _data->get_channels()) {
181 size_t n_gaussians_c = 0;
183 auto& gaussians_src = gaussians_srcs[channel];
186 n_gaussians_c += gaussians_src.
back()->size();
188 if (n_gaussians_src == 0) {
189 n_gaussians_src = n_gaussians_c;
190 }
else if (n_gaussians_src != n_gaussians_c) {
195 return gaussians_srcs;
205 double _evaluate_observation(
size_t idx,
bool print =
false,
bool skip_zeroing =
false) {
208 const auto& channel = _data->at(idx).
get().get_channel();
209 const size_t n_gaussians_psf = _psfcomps[idx]->
size();
210 const auto& gaussians_src = _gaussians_srcs.
at(channel);
212 const auto& factors_extra_in = _factors_extra_in[idx].lock();
213 const auto& factors_grad_in = _factors_grad_in[idx].lock();
216 for (
size_t i_psf = 0; i_psf < n_gaussians_psf; ++i_psf) {
220 std::cout <<
"Setting factors for psf[" << i_psf <<
"], src[" << i_src <<
"]"
222 const size_t n_gauss_src = gaussians_src[i_src++]->
size();
223 if (n_gauss_src !=
src->get_n_gaussians(channel)) {
225 this->
str() +
"._data[" + channel.str() +
"].get_n_gaussians(channel)="
227 +
" != get_gaussians(channel).size()=" +
std::to_string(n_gauss_src));
229 src->set_grad_param_factors(channel, _factors_grad, offset);
230 src->set_extra_param_factors(channel, _factors_extra, offset);
232 const size_t row_max =
offset + n_gauss_src;
233 for (
size_t row = offset;
row < row_max; ++
row) {
238 factors_extra_in->set_value_unchecked(
row,
col, _factors_extra[
row][
col]);
247 factors_grad_in->set_value_unchecked(
row,
col, _factors_grad[
row][
col]);
260 auto& grads = *(this->_grads[idx]);
261 for (
auto& grad : grads) grad->
fill(0);
264 return _evaluators[idx]->loglike_pixel();
273 double _evaluate_priors(
bool print =
false,
bool normalize_loglike =
false) {
276 size_t idx_resid = 0;
279 for (
const auto& prior : _priors) {
280 auto eval = prior->evaluate(is_jacobian, normalize_loglike);
282 if (_residuals_prior !=
nullptr) {
283 size_t idx_jac = idx_resid;
284 size_t n_resid = eval.residuals.size();
285 for (
const double value : eval.residuals) {
286 this->_residuals_prior->set_value(0, idx_resid++, value);
291 for (
const auto& [param, values] : eval.jacobians) {
292 if (values.size() != n_resid) {
293 throw std::logic_error(
"Prior=" + prior->str() +
" param=" + param.get().str()
297 const size_t idx_param = _offsets_params.at(param);
298 const auto& img = _outputs_prior.
at(idx_param);
299 for (
const double value : values) {
300 img->set_value(0, idx_jac++, value);
305 idx_resid += n_resid;
308 std::cout <<
"Evaluated prior with loglike=" << eval.loglike <<
"; at idx_resid=" << idx_resid
332 _check_obs_idx(idx_obs);
335 const Channel& channel = observation.get_channel();
336 const auto& psfcomps = _psfcomps.
at(idx_obs);
337 const size_t n_c = observation.get_n_cols();
338 const size_t n_r = observation.get_n_rows();
340 const size_t n_outputs =
outputs.size();
341 const bool has_outputs = n_outputs > 0;
342 const auto* outputs_obs = has_outputs ? &(
outputs.at(idx_obs)) : nullptr;
350 if (is_mode_jacobian || is_mode_loglike_grad) {
351 const size_t size_out = outputs_obs->size();
352 size_t n_obsect = is_mode_jacobian ? (_n_params_free + 1) : this->
size();
353 if (size_out != n_obsect) {
356 + (is_mode_jacobian ?
"free+1" :
"obs") +
"="
360 if (!(n_outputs > idx_obs)) {
368 = !do_output ? nullptr
369 : (has_outputs ? outputs_obs->at(0)
370 : std::make_shared<Image>(n_r, n_c,
nullptr, coordsys));
372 auto data_eval = !is_mode_image ? observation.get_image_ptr_const() :
nullptr;
373 auto sigma_inv = !is_mode_image ? observation.get_sigma_inverse_ptr_const() :
nullptr;
375 const auto& gaussians_src = _gaussians_srcs.
at(channel);
376 const size_t n_gaussians_psf = psfcomps->size();
378 for (
size_t p = 0; p < n_gaussians_psf; ++p) {
379 const auto& psfcomp = psfcomps->at(p);
380 const auto value = psfcomp.get_integral_value();
388 size_t n_gaussians_src = 0;
389 for (
const auto& gaussians : gaussians_src) n_gaussians_src +=
gaussians->
size();
390 const size_t n_gaussians_conv = n_gaussians_src * n_gaussians_psf;
394 for (
auto psfit = psfcomps->cbegin(); psfit < psfcomps->cend(); psfit++) {
395 for (
const auto& gaussians : gaussians_src) {
396 for (
const auto& gaussian : *
gaussians) {
401 auto gaussians = std::make_shared<const ConvolvedGaussians>(
data);
402 if (
gaussians->size() != n_gaussians_conv) {
413 if (is_mode_jacobian || is_mode_loglike_grad) {
420 this->_make_param_maps<true>(factors_extra_weak, factors_grad_weak, map_extra_weak,
421 map_grad_weak, n_gaussians_conv, n_gaussians_psf, channel,
422 coordsys, idx_obs, map_extra_in, map_grad_in, factors_extra_in,
425 this->_make_param_maps<false>(factors_extra_weak, factors_grad_weak, map_extra_weak,
426 map_grad_weak, n_gaussians_conv, n_gaussians_psf, channel,
427 coordsys, idx_obs, map_extra_in, map_grad_in, factors_extra_in,
433 if (is_mode_jacobian) {
435 const size_t n_free = _n_params_free + 1;
437 for (
size_t idx_img = 0; idx_img < n_free; ++idx_img) {
440 image = (*outputs_obs)[idx_img];
441 if (
image ==
nullptr) {
445 const auto& image_ref = observation.get_image();
451 +
" incompatible with corresponding image="
453 +
" due to: " + msg);
456 image = std::make_shared<Image>(n_r, n_c,
nullptr, coordsys);
460 grads = std::make_shared<ImageArray<T, Image>>(&
arrays);
461 }
else if (is_mode_loglike_grad) {
465 <<
"] to image of n_cols=_n_params_free + 1=" << _n_params_free + 1
466 <<
"; _offsets_params.size()=" << _offsets_params.size() <<
std::endl;
468 auto image = has_outputs ? outputs_obs->at(0)
469 : std::make_shared<Image>(1, _n_params_free + 1,
nullptr, coordsys);
471 grads = std::make_shared<ImageArray<T, Image>>(&
arrays);
478 std::cout <<
"output is null: " << (output ==
nullptr) <<
"\n";
483 auto output2 = output;
487 auto test =
Evaluator(gaussians, data_eval, sigma_inv, output, residual, grads, map_grad_in,
488 factors_grad_in, map_extra_in, factors_extra_in);
490 = std::make_unique<Evaluator>(gaussians, data_eval, sigma_inv, output, residual, grads,
491 map_grad_in, factors_grad_in, map_extra_in, factors_extra_in
498 template <
bool pr
int>
501 size_t n_gaussians_conv,
size_t n_gaussians_psf,
const Channel& channel,
507 unsigned int n_obsired = map_extra_weak.
expired() + map_grad_weak.
expired()
509 const bool expired = n_obsired == 4;
511 std::cout <<
"expired=" << expired <<
"\n";
513 if (!((n_obsired == 0) || expired)) {
516 auto map_extra_mut = expired ? std::make_shared<Indices>(n_gaussians_conv, 2,
nullptr, coordsys)
517 : map_extra_weak.
lock();
518 auto map_grad_mut = expired ? std::make_shared<Indices>(
520 : map_grad_weak.
lock();
521 auto factors_extra_mut = expired ? std::make_shared<Image>(n_gaussians_conv, 3,
nullptr, coordsys)
522 : factors_extra_weak.
lock();
523 auto factors_grad_mut = expired ? std::make_shared<Image>(
525 : factors_grad_weak.
lock();
532 factors_extra.
clear();
533 factors_grad.clear();
535 map_extra.reserve(n_gaussians_conv);
536 map_grad.reserve(n_gaussians_conv);
537 factors_extra.reserve(n_gaussians_conv);
538 factors_grad.reserve(n_gaussians_conv);
545 for (
size_t i_psf = 0; i_psf < n_gaussians_psf; ++i_psf) {
547 src->add_grad_param_map(channel, map_grad, _offsets_params);
548 src->add_extra_param_map(channel, map_extra, map_grad, _offsets_params);
550 src->add_grad_param_factors(channel, factors_grad);
551 src->add_extra_param_factors(channel, factors_extra);
553 if constexpr (print) {
555 std::cout << str_map_ref<true>(_offsets_params) <<
std::endl;
558 if (map_grad.size() != n_gaussians_conv) {
577 if (sizes != sizes_mut) {
579 "make_evaluator Jacobian map_extra,map_grad,factors_extra,factors_grad sizes_new="
580 + sizes +
"!=sizes_mut=" + sizes_mut +
"; did you make a new evaluator with different"
581 "free parameters from an old one?");
588 for (
size_t idx_row = 0; idx_row < map_extra.size(); idx_row++) {
589 const auto&
row = map_extra[idx_row];
591 const auto value_mut = map_extra_mut->get_value_unchecked(idx_row,
col);
592 if ((idx_err < 2) && (
row[
col] != value_mut))
593 errmsgs[idx_err++] =
"map_extra[" +
std::to_string(idx_row) +
"]["
597 if constexpr (print) {
603 if constexpr (print) {
606 for (
size_t idx_row = 0; idx_row < map_grad.size(); idx_row++) {
607 const auto&
row = map_grad[idx_row];
609 const auto value_mut = map_grad_mut->get_value_unchecked(idx_row,
col);
610 if ((idx_err < 4) && (
row[
col] != value_mut))
615 if constexpr (print) {
622 std::string errmsg_extra = errmsgs[0] + errmsgs[1];
624 if (!errmsg_extra.
empty() || !errmsg_grad.
empty()) {
626 =
"make_evaluator Jacobian map_extra/map_grad have different values"
627 " from old; did you make a new evaluator with different free parameters from an "
628 " old one? Errors:\n"
630 if (!errmsg_extra.
empty() && !errmsg_grad.
empty()) errmsg +=
"...\n";
631 errmsg += errmsg_grad;
639 if constexpr (print) {
642 for (
size_t row = 0;
row < n_gaussians_conv; ++
row) {
644 map_extra_mut->set_value_unchecked(
row,
col, (map_extra)[
row][
col]);
647 factors_extra_mut->set_value_unchecked(
row,
col, (factors_extra)[
row][
col]);
650 map_grad_mut->set_value_unchecked(
row,
col, (map_grad)[
row][
col]);
651 factors_grad_mut->set_value_unchecked(
row,
col, (factors_grad)[
row][
col]);
653 if constexpr (print) {
667 _factors_extra_in[idx_obs] = factors_extra_mut;
668 _factors_grad_in[idx_obs] = factors_grad_mut;
672 factors_extra_in =
std::move(factors_extra_mut);
673 factors_grad_in =
std::move(factors_grad_mut);
679 for (
auto& source : _sources) source->add_extra_param_map(channel, map_extra, map_grad, offsets);
682 for (
auto& source : _sources) source->add_extra_param_factors(channel, factors);
685 for (
auto& source : _sources) source->add_grad_param_map(channel, map, offsets);
688 for (
auto& source : _sources) source->add_grad_param_factors(channel, factors);
705 bool verify =
false,
double findiff_frac = 1e-5,
706 double findiff_add = 1e-5,
double rtol = 1e-3,
707 double atol = 1e-8) {
713 params_free = nonconsecutive_unique<ParamBaseCRef>(params_free);
715 const size_t n_params = this->_offsets_params.size();
716 if (n_params != params_free.size()) {
719 +
" != params_free=" + str_iter_ref<true>(params_free) +
".size() + 1 ="
723 for (
size_t idx_param = 0; idx_param < n_params; idx_param++) {
724 const auto param = params_free[idx_param];
725 if (_offsets_params.find(param) == _offsets_params.end()) {
729 param_idx[param] = idx_param;
739 for (
const auto& grads_obs : _grads) {
740 const auto& values = grads_obs->at(0);
741 const size_t n_cols = values.get_n_cols();
742 if (n_cols != n_params + 1) {
745 +
" != n_params (from offsets) = " +
std::to_string(n_params + 1));
747 for (
const auto& [paramref,
col] : _offsets_params) {
748 loglike_grads.
at(param_idx[paramref]) += values.get_value(0,
col);
753 for (
const auto& prior : _priors) {
754 auto result = prior->evaluate(
true);
757 for (
const auto& paramref : params_free) {
758 double dll_dx =
result.compute_dloglike_dx(paramref);
759 loglike_grads.
at(param_idx[paramref]) += dll_dx;
769 params = nonconsecutive_unique<ParamBaseRef>(params);
771 for (
size_t idx_param = 0; idx_param < n_params; idx_param++) {
772 auto& param = params[idx_param].get();
773 double value = param.get_value_transformed();
774 double diff =
std::copysign(std::abs(value * findiff_frac) + findiff_add,
775 loglike_grads[idx_param]);
780 param.set_value_transformed(value + diff / 2.);
781 loglike_new_plus = this->
evaluate();
782 param.set_value_transformed(value - diff / 2.);
783 loglike_new_minus = this->
evaluate();
785 param.set_value_transformed(value + diff);
786 loglike_new_plus = this->
evaluate();
789 auto dloglike = (
sum_iter(loglike_new_plus) -
sum_iter(loglike_new_minus)) / diff;
790 auto close =
isclose(dloglike, loglike_grads[idx_param], rtol, atol);
791 param.set_value_transformed(value);
792 if (!close.isclose) {
793 double dll_dx_exact = loglike_grads[idx_param];
794 double dlp_dx_sum = 0;
795 for (
const auto& prior : _priors) {
796 double dlp_dx = prior->evaluate(
true).compute_dloglike_dx(param,
true);
797 dlp_dx_sum += dlp_dx;
799 double dlp_dx_findiff = (loglike_new_plus.
back() - loglike_new_minus.
back()) / diff;
800 errmsg += param.str() +
" failed loglike_grad verification; isclose=" + close.str()
812 return loglike_grads;
826 std::optional<HessianOptions>
options = std::nullopt,
827 bool print =
false) {
832 params_free = nonconsecutive_unique<ParamBaseRef>(params_free);
833 const size_t n_params_free = params_free.size();
836 transforms.
reserve(n_params_free);
838 for (
auto& paramref : params_free) {
839 auto& param = paramref.get();
840 transforms.emplace_back(param.get_transform_ptr());
841 double value = param.get_value();
842 param.set_transform(
nullptr);
843 if (param.get_value() != value) {
850 auto hessian = std::make_unique<Image>(n_params_free, n_params_free);
851 const auto n_obs = this->
size();
854 size_t n_integral = 0;
857 for (
size_t idx_param = 0; idx_param < n_params_free; idx_param++) {
858 const auto param = params_free[idx_param];
859 const auto found = _offsets_params.find(param);
860 if (found == _offsets_params.end()) {
864 param_idx[idx_param] = (*found).second;
866 const auto name = param.get().get_name();
877 size_t idx_param = 0;
878 for (
auto& paramref : params_free) {
879 auto& param = paramref.get();
880 const double value_init = param.get_value();
881 const double value_transformed = param.get_value_transformed();
887 double diffabs = std::abs(diff);
891 for (
size_t idx_col = 0; idx_col < n_params_free; ++idx_col) {
892 double dll_dx = loglike_grad2[idx_col] -
loglike_grad[idx_col];
893 dll_dx *= opt.
return_negative ? (1 - 2 * (dll_dx > 0)) / diffabs : 1 / diff;
894 hessian->set_value_unchecked(idx_param, idx_col, dll_dx);
896 param.set_value(value_init);
900 if ((n_integral > 0) && (n_frac > 0)) {
902 "Cannot compute_hessian without options for model with free Integral and "
904 "parameters (this is a bug); please supply options to use finite differencing.");
912 for (
size_t idx_obs = 0; idx_obs < n_obs; ++idx_obs) {
913 const auto& jacs = this->_grads.
at(idx_obs);
914 for (
size_t idx_param1 = 0; idx_param1 < n_params_free; ++idx_param1) {
915 const auto& jac1 = jacs->at(param_idx[idx_param1]);
916 const size_t n_rows = jac1.get_n_rows();
917 const size_t n_cols = jac1.get_n_cols();
919 for (
size_t idx_param2 = 0; idx_param2 < n_params_free; ++idx_param2) {
920 const auto& jac2 = jacs->at(param_idx[idx_param2]);
921 const size_t n_rows2 = jac2.get_n_rows();
922 const size_t n_cols2 = jac2.get_n_cols();
924 if ((n_rows != n_rows2) || (n_cols != n_cols2)) {
931 for (
size_t row = 0;
row < n_rows; ++
row) {
932 for (
size_t col = 0;
col < n_cols; ++
col) {
935 dx -= jac1.get_value(
row,
col) * jac2.get_value(
row,
col);
939 std::cout <<
"calling hessian->add_value(" << idx_param1 <<
"," << idx_param2
942 hessian->add_value(idx_param1, idx_param2, dx);
946 if (include_prior && (_priors.
size() > 0)) {
948 for (
size_t idx_param = 0; idx_param < n_params_free; idx_param++) {
949 const auto param = params_free[idx_param];
950 if (_offsets_params.find(param) == _offsets_params.end()) {
954 param_idx_map[param] = idx_param;
957 size_t prior_size = 0;
958 for (
const auto& prior : _priors) {
959 prior_size += prior->
size();
961 auto dll_prior = std::make_unique<Image>(n_params_free, prior_size);
964 size_t idx_prior = 0;
965 for (
const auto& prior : _priors) {
966 auto result = prior->evaluate(
true);
968 for (
const auto& [param_cref, values_base] :
result.jacobians) {
969 size_t idx_residual = 0;
970 size_t idx_param = param_idx_map[param_cref];
971 for (
const auto& value_base : values_base) {
972 dll_prior->add_value(idx_param, idx_prior + idx_residual++, value_base);
975 idx_prior += prior->
size();
978 for (
size_t row = 0;
row < n_params_free; ++
row) {
979 for (
size_t row2 = 0; row2 < n_params_free; ++row2) {
981 for (
size_t col = 0;
col < prior_size; ++
col) {
982 dll_dx -= dll_prior->get_value(
row,
col) * dll_prior->get_value(row2,
col);
984 hessian->add_value(
row, row2, dll_dx);
986 std::cout <<
"calling hessian->add_value(" <<
row <<
"," << row2 <<
"," << dll_dx
995 size_t idx_param = 0;
996 for (
auto& paramref : params_free) {
997 auto& param = paramref.get();
998 double value = param.get_value();
999 param.set_transform(transforms[idx_param]);
1000 if (param.get_value() != value) {
1020 if (!_is_setup)
throw std::runtime_error(
"Can't call evaluate before setup_evaluators");
1025 if (is_loglike_grad) {
1029 for (
size_t idx = 0; idx < _size; ++idx) {
1030 this->_grads[idx]->
at(0).fill(0);
1034 for (
size_t idx = 0; idx < _size; ++idx) {
1035 result[idx] = this->_evaluate_observation(idx, print, is_loglike_grad);
1039 result[_size] = this->_evaluate_priors(print, normalize_loglike);
1052 _check_obs_idx(idx);
1053 return _evaluate_observation(_evaluators[idx]);
1066 for (
auto& source : _sources) {
1067 auto data = source->get_gaussians(channel)->get_data();
1068 in.emplace_back(
data);
1071 return std::make_unique<lsst::gauss2d::Gaussians>(in);
1080 const size_t n_data = this->
size();
1081 if (this->_likelihood_const_terms.
empty()) {
1082 this->_likelihood_const_terms.
resize(n_data + 1);
1084 for (
size_t idx = 0; idx < n_data; ++idx) {
1085 const auto& sigma_inv = this->_data[idx].get_sigma_inverse();
1087 const size_t n_rows = sigma_inv.get_n_rows();
1088 const size_t n_cols = sigma_inv.get_n_cols();
1089 for (
size_t col = 0;
col < n_cols; ++
col) {
1090 for (
size_t row = 0;
row < n_rows; ++
row) {
1094 _likelihood_const_terms[idx] =
loglike;
1098 for (
const auto& prior : this->_priors) {
1099 auto values = prior->get_loglike_const_terms();
1100 for (
const auto value : values)
loglike += value;
1102 _likelihood_const_terms[n_data] =
loglike;
1103 return _likelihood_const_terms;
1108 for (
auto& source : _sources) n_g += source->get_n_gaussians(channel);
1117 for (
auto& [param, offset] : _offsets_params) {
1118 offsets.emplace_back(param, offset);
1124 for (
auto& source : _sources) source->get_parameters(params, filter);
1125 for (
auto& psfmodel : _psfmodels) psfmodel->get_parameters(params, filter);
1126 _data->get_parameters(params, filter);
1131 for (
auto& source : _sources) source->get_parameters_const(params, filter);
1132 for (
auto& psfmodel : _psfmodels) psfmodel->get_parameters_const(params, filter);
1133 _data->get_parameters_const(params, filter);
1140 _check_obs_idx(idx);
1141 for (
auto& source : _sources) source->get_parameters(params, filter);
1142 _psfmodels[idx]->get_parameters(params, filter);
1143 _data->at(idx).
get().get_parameters(params, filter);
1150 _check_obs_idx(idx);
1151 for (
auto& source : _sources) source->get_parameters_const(params, filter);
1152 _psfmodels[idx]->get_parameters_const(params, filter);
1153 _data->at(idx).
get().get_parameters_const(params, filter);
1167 size_t index)
const override {
1168 for (
auto& source : _sources) {
1169 source->set_extra_param_factors(channel, factors, index);
1170 index += source->get_n_gaussians(channel);
1175 size_t index)
const override {
1176 for (
auto& source : _sources) {
1177 source->set_grad_param_factors(channel, factors, index);
1178 index += source->get_n_gaussians(channel);
1204 bool print =
false) {
1205 const size_t n_outputs = outputs.size();
1206 const bool has_outputs = n_outputs > 0;
1208 if (n_outputs !=
size()) {
1213 const size_t n_residuals = residuals.
size();
1214 const bool has_residuals = n_residuals > 0;
1215 if (has_residuals) {
1216 if (n_residuals !=
size()) {
1221 if (force || (!_is_setup || (mode != _mode))) {
1222 _evaluators.
clear();
1225 _outputs_prior.
clear();
1226 _residuals_prior =
nullptr;
1234 if (is_jacobian || is_loglike_grad) {
1237 params_free = nonconsecutive_unique<ParamBaseCRef>(params_free);
1238 _n_params_free = params_free.
size();
1240 std::cout <<
"n_params_free=" << _n_params_free <<
" from params_free:" <<
std::endl;
1241 std::cout << str_iter_ref<true>(params_free) <<
std::endl;
1245 _offsets_params.clear();
1248 size_t n_residuals_prior = 0;
1249 bool do_residuals_prior = (_size_priors > 0)
1252 if (do_residuals_prior) {
1253 for (
const auto& prior : this->_priors) {
1254 n_residuals_prior += prior->size();
1258 = residuals_prior ==
nullptr
1259 ? (_residuals_prior = std::make_shared<Image>(1, n_residuals_prior))
1261 if ((_residuals_prior->get_n_rows() != 1)
1262 || (_residuals_prior->get_n_cols() != n_residuals_prior)) {
1267 for (
size_t col = 0;
col < n_residuals_prior; ++
col) {
1268 _residuals_prior->set_value(0,
col, 0);
1278 const size_t n_params_jac = _n_params_free + 1;
1279 if (outputs_prior.
size() == 0) {
1280 for (
size_t idx_jac = 0; idx_jac < n_params_jac; ++idx_jac) {
1281 _outputs_prior.
emplace_back(std::make_shared<Image>(1, n_residuals_prior));
1284 }
else if (outputs_prior.
size() != n_params_jac) {
1287 +
" != (n_params_free + 1 = " +
std::to_string(n_params_jac) +
")");
1289 _outputs_prior.
reserve(n_params_jac);
1291 for (
auto& output_prior : outputs_prior) {
1292 if (output_prior ==
nullptr) {
1296 if ((output_prior->get_n_rows() != 1)
1297 || (output_prior->get_n_cols() != n_residuals_prior)) {
1299 +
"]=" + output_prior->str() +
" rows,cols != 1,"
1302 for (
size_t col = 0;
col < n_residuals_prior; ++
col) {
1303 output_prior->set_value(0,
col, 0);
1315 const auto&
channels = _data->get_channels();
1319 for (
size_t idx = 0; idx < _size; ++idx) {
1320 auto result = this->_make_evaluator(idx, mode, outputs,
1321 has_residuals ? residuals[idx] : nullptr, print);
1323 if (is_jacobian || is_loglike_grad) {
1333 size_t idx_prior = 0;
1336 for (
const auto& prior : _priors) {
1337 auto eval_prior = prior->evaluate(is_jacobian);
1338 size_t n_resid = prior->
size();
1339 size_t n_resid_eval = eval_prior.residuals.size();
1340 if (n_resid != eval_prior.residuals.size()) {
1342 +
" != evaluate(true).residuals.size()="
1346 size_t idx_param = 0;
1347 for (
const auto& param_jacob : eval_prior.jacobians) {
1348 const auto& param = param_jacob.first;
1349 if (params_prior.
find(param) != params_prior.
end()) {
1353 +
"]=" + param.get().str() +
" already in prior param set."
1354 +
" Did you apply multiple priors to the same parameter?");
1362 if (_size_priors > 0) {
1363 if ((_residuals_prior !=
nullptr) && (_residuals_prior->get_n_cols() != n_residuals_prior)) {
1375 size_t size()
const {
return _size; }
1379 std::string str = type_name_str<Model>(
false, namespace_separator) +
"("
1380 + (name_keywords ?
"data=" :
"") + _data->repr(name_keywords, namespace_separator)
1381 +
", " + (name_keywords ?
"psfmodels=[" :
"[");
1382 for (
const auto&
x : _psfmodels) str +=
x->repr(name_keywords, namespace_separator) +
",";
1383 str += (name_keywords ?
"], sources=[" :
"], [");
1384 for (
const auto& s : _sources) str += s->repr(name_keywords, namespace_separator) +
",";
1385 str += (name_keywords ?
"], priors=[" :
"], [");
1390 std::string str = type_name_str<Model>(
true) +
"(data=" + _data->str() +
", psfmodels=[";
1391 for (
const auto&
x : _psfmodels) str +=
x->str() +
",";
1392 str +=
"], sources=[";
1393 for (
const auto&
x : _sources) str +=
x->str() +
",";
1394 str +=
"], priors=[";
1395 for (
const auto&
x : _priors) str +=
x->str() +
",";
1413 double rtol = 1e-3,
double atol = 1e-3,
double max_ll_diff = 0) {
1414 if (!(max_ll_diff >= 0)) {
1421 const size_t n_obs = this->
size();
1425 for (
size_t idx = 0; idx < n_obs; ++idx) {
1426 const auto& grads = *(_grads.
at(idx));
1427 const auto& observation = _data->at(idx).
get();
1428 const auto& sigma_inv = observation.get_sigma_inverse();
1430 const Channel& channel = observation.get_channel();
1431 const size_t n_cols = observation.get_n_cols();
1432 const size_t n_rows = observation.get_n_rows();
1434 filter.channel = channel;
1440 size_t idx_param_max = 0;
1441 for (
const auto& param : params) {
1442 auto found = _offsets_params.find(param);
1443 if (found == _offsets_params.end()) {
1446 +
" not found in offsets; was it freed after setup_evaluators was called?");
1448 size_t idx_param = found->second;
1449 if (idx_param > idx_param_max) idx_param_max = idx_param;
1452 if (!(idx_param_max < grads.size())) {
1455 +
"; is the jacobian array large enough?");
1458 double loglike_jac = this->_evaluate_observation(idx);
1459 const auto& grad = grads[0];
1460 size_t n_nonzero = 0;
1461 for (
unsigned int i = 0; i < n_cols; ++i) {
1462 for (
unsigned int j = 0; j < n_rows; ++j) {
1463 n_nonzero += grad.get_value(j, i) != 0;
1466 if (n_nonzero > 0) {
1471 double loglike_img =
result.first->loglike_pixel();
1472 auto is_close_ll =
isclose(loglike_img, loglike_jac, 0., max_ll_diff);
1473 if (!is_close_ll.isclose) {
1475 +
to_string_float(loglike_img) +
" !close (isclose=" + is_close_ll.str()
1479 const auto& image_current = *(
result.second.first);
1481 for (
auto& param_ref : params) {
1482 auto& param = param_ref.get();
1483 const size_t idx_param = _offsets_params.find(param_ref)->second;
1484 const auto& grad = grads[idx_param];
1486 const double value_init = param.get_value();
1487 const double value_transformed = param.get_value_transformed();
1488 double diff = value_transformed * findiff_frac;
1489 if (std::abs(diff) < findiff_add) diff = findiff_add;
1494 result_diff.first->loglike_pixel();
1496 size_t n_failed = 0;
1498 for (
unsigned int i = 0; i < n_cols; ++i) {
1499 for (
unsigned int j = 0; j < n_rows; ++j) {
1500 double delta = sigma_inv.get_value(j, i) / diff
1501 * (result_diff.second.first->get_value(j, i)
1502 - image_current.get_value(j, i));
1503 double jac = grad.get_value(j, i);
1505 auto close =
isclose(jac, delta, rtol, atol);
1506 if (!close.isclose) {
1513 param.set_value(value_init);
1514 const double value_new = param.get_value();
1515 if (value_new != value_init) {
1522 double median = ratios[ratios.
size() / 2];
1530 filter.channel = std::nullopt;
1536 for (
const auto& prior : this->_priors) {
1537 auto result_base = prior->evaluate(
true);
1538 auto residuals_base = result_base.residuals;
1539 const size_t n_values = residuals_base.size();
1541 for (
const auto& [param_cref, values_base] : result_base.jacobians) {
1542 if (values_base.size() != n_values) {
1549 auto found =
std::find(params.begin(), params.end(), param_cref);
1550 if (found == params.end()) {
1551 throw std::logic_error(
"Couldn't find " + param_cref.get().str() +
" in all params");
1553 auto& param = (*found).get();
1555 const double value_init = param.get_value();
1556 const double value_transformed = param.get_value_transformed();
1557 double delta = value_transformed * findiff_frac;
1558 if (std::abs(delta) < findiff_add) delta = findiff_add;
1560 auto result_new = prior->evaluate(
true);
1562 const size_t idx_param = _offsets_params.at(param);
1565 size_t n_failed = 0;
1566 size_t idx_value = 0;
1567 for (
const double jac_base : values_base) {
1569 = (result_new.residuals[idx_value] - result_base.residuals[idx_value]) / delta;
1570 auto close =
isclose(jac_base, jac_findiff, rtol, atol);
1571 if (!close.isclose) {
1573 prior->evaluate(
true);
1574 ratios.
push_back(jac_base / jac_findiff);
1578 param.set_value(value_init);
1579 const double value_new = param.get_value();
1580 if (value_new != value_init) {
1587 double median = ratios[ratios.
size() / 2];
1589 +
" failed for prior=" + prior->str()