SourceXtractorPlusPlus  0.16
Please provide a description of the project.
ModelFittingConfig.cpp
Go to the documentation of this file.
1 
17 /*
18  * @file ModelFittingConfig.cpp
19  * @author Nikolaos Apostolakos <nikoapos@gmail.com>
20  */
21 
22 #include <string>
23 #include <boost/python/extract.hpp>
24 #include <boost/python/object.hpp>
25 #include <boost/python/tuple.hpp>
26 #include <boost/python/dict.hpp>
27 
28 #include "ElementsKernel/Logging.h"
29 
31 
37 #include "Pyston/GIL.h"
38 #include "Pyston/Exceptions.h"
41 
42 #include <string>
43 #include <boost/python/extract.hpp>
44 #include <boost/python/object.hpp>
45 
46 #ifdef WITH_ONNX_MODELS
48 #endif
49 
51 
52 
53 using namespace Euclid::Configuration;
54 namespace py = boost::python;
56 
58 
59 namespace SourceXtractor {
60 
61 template<typename Signature>
63 };
64 
65 template<>
66 struct FunctionFromPython<double(const SourceInterface&)> {
67  static
68  std::function<double(const SourceInterface&)> get(const char *readable,
70  py::object py_func) {
71  auto wrapped = builder.build<double(const AttributeSet&)>(py_func, ObjectInfo{});
72 
73  if (!wrapped.isCompiled()) {
74  logger.warn() << "Could not compile " << readable << ": " << wrapped.reason()->what();
75  wrapped.reason()->log(log4cpp::Priority::DEBUG, logger);
76  }
77  else {
78  logger.info() << readable << " compiled";
79  Pyston::GraphvizGenerator gv(readable);
80  wrapped.getTree()->visit(gv);
81  logger.info() << gv.str();
82  }
83 
84  return [wrapped](const SourceInterface& o) -> double {
85  return wrapped(ObjectInfo{o});
86  };
87  }
88 };
89 
90 template<>
91 struct FunctionFromPython<double(const Pyston::Context&, const std::vector<double>&)> {
92  static
93  std::function<double(const Pyston::Context&, const std::vector<double>&)>
94  get(const char *readable, Pyston::ExpressionTreeBuilder& builder, py::object py_func,
95  size_t nparams) {
96  auto wrapped = builder.build<double(const std::vector<double>&)>(py_func, nparams);
97  if (!wrapped.isCompiled()) {
98  logger.warn() << "Could not compile " << readable << ": " << wrapped.reason()->what();
99  wrapped.reason()->log(log4cpp::Priority::DEBUG, logger);
100  }
101  else {
102  logger.info() << readable << " compiled";
103  Pyston::GraphvizGenerator gv(readable);
104  wrapped.getTree()->visit(gv);
105  logger.info() << gv.str();
106  }
107 
108  return wrapped;
109  }
110 };
111 
112 template<>
113 struct FunctionFromPython<double(double, const SourceInterface&)> {
114  static
115  std::function<double(double, const SourceInterface&)> get(const char *readable,
117  py::object py_func) {
118  auto wrapped = builder.build<double(double, const AttributeSet&)>(py_func, ObjectInfo{});
119 
120  if (!wrapped.isCompiled()) {
121  logger.warn() << "Could not compile " << readable << ": " << wrapped.reason()->what();
122  wrapped.reason()->log(log4cpp::Priority::DEBUG, logger);
123  }
124  else {
125  logger.info() << readable << " compiled";
126  Pyston::GraphvizGenerator gv(readable);
127  wrapped.getTree()->visit(gv);
128  logger.info() << gv.str();
129  }
130 
131  return [wrapped](double a, const SourceInterface& o) -> double {
132  return wrapped(a, ObjectInfo{o});
133  };
134  }
135 };
136 
137 ModelFittingConfig::ModelFittingConfig(long manager_id) : Configuration(manager_id) {
138  declareDependency<PythonConfig>();
139 }
140 
142  Pyston::GILLocker locker;
143  m_parameters.clear();
144  m_models.clear();
145  m_frames.clear();
146  m_priors.clear();
147  m_outputs.clear();
148 }
149 
151  Pyston::GILLocker locker;
152  try {
153  initializeInner();
154  }
155  catch (py::error_already_set &e) {
156  throw Pyston::Exception().log(log4cpp::Priority::ERROR, logger);
157  }
158 }
159 
160 static double image_to_world_alpha(const Pyston::Context& context, double x, double y) {
161  auto coord_system = boost::any_cast<std::shared_ptr<CoordinateSystem>>(context.at("coordinate_system"));
162  return coord_system->imageToWorld({x, y}).m_alpha;
163 }
164 
165 static double image_to_world_delta(const Pyston::Context& context, double x, double y) {
166  auto coord_system = boost::any_cast<std::shared_ptr<CoordinateSystem>>(context.at("coordinate_system"));
167  return coord_system->imageToWorld({x, y}).m_delta;
168 }
169 
171  Pyston::ExpressionTreeBuilder expr_builder;
172  expr_builder.registerFunction<double(const Pyston::Context&, double, double)>(
173  "image_to_world_alpha", image_to_world_alpha
174  );
175  expr_builder.registerFunction<double(const Pyston::Context&, double, double)>(
176  "image_to_world_delta", image_to_world_delta
177  );
178 
179  /* Constant parameters */
180  for (auto& p : getDependency<PythonConfig>().getInterpreter().getConstantParameters()) {
181  auto value_func = FunctionFromPython<double(const SourceInterface&)>::get(
182  "Constant parameter", expr_builder, p.second.attr("get_value")
183  );
184 
185  m_parameters[p.first] = std::make_shared<FlexibleModelFittingConstantParameter>(
186  p.first, value_func);
187  }
188 
189  /* Free parameters */
190  for (auto& p : getDependency<PythonConfig>().getInterpreter().getFreeParameters()) {
191  auto init_value_func = FunctionFromPython<double(const SourceInterface&)>::get(
192  "Free parameter", expr_builder, p.second.attr("get_init_value")
193  );
194 
195  auto py_range_obj = p.second.attr("get_range")();
196 
198  std::string type_string(py::extract<char const*>(py_range_obj.attr("__class__").attr("__name__")));
199 
200  if (type_string == "Unbounded") {
201  auto factor_func = FunctionFromPython<double(double, const SourceInterface&)>::get(
202  "Unbounded", expr_builder, py_range_obj.attr("get_normalization_factor")
203  );
204  converter = std::make_shared<FlexibleModelFittingUnboundedConverterFactory>(factor_func);
205  } else if (type_string == "Range") {
206  auto min_func = FunctionFromPython<double(double, const SourceInterface&)>::get(
207  "Range min", expr_builder, py_range_obj.attr("get_min")
208  );
209  auto max_func = FunctionFromPython<double(double, const SourceInterface&)>::get(
210  "Range max", expr_builder, py_range_obj.attr("get_max")
211  );
212 
213  auto range_func = [min_func, max_func] (double init, const SourceInterface& o) -> std::pair<double, double> {
214  return {min_func(init, o), max_func(init, o)};
215  };
216 
217  bool is_exponential = py::extract<int>(py_range_obj.attr("get_type")().attr("value")) == 2;
218 
219  if (is_exponential) {
220  converter = std::make_shared<FlexibleModelFittingExponentialRangeConverterFactory>(range_func);
221  } else {
222  converter = std::make_shared<FlexibleModelFittingLinearRangeConverterFactory>(range_func);
223  }
224  } else {
225  throw Elements::Exception("Unknown converter type: " + type_string);
226  }
227  m_parameters[p.first] = std::make_shared<FlexibleModelFittingFreeParameter>(
228  p.first, init_value_func, converter);
229  }
230 
231  /* Dependent parameters */
232  for (auto& p : getDependency<PythonConfig>().getInterpreter().getDependentParameters()) {
233  auto py_func = p.second.attr("func");
235  py::list param_ids = py::extract<py::list>(p.second.attr("params"));
236  for (int i = 0; i < py::len(param_ids); ++i) {
237  int id = py::extract<int>(param_ids[i]);
238  params.push_back(m_parameters[id]);
239  }
240 
241  auto dependent = FunctionFromPython<double(const Pyston::Context&, const std::vector<double>&)>
242  ::get("Dependent parameter", expr_builder, py_func, params.size());
243 
244  auto dependent_func = [dependent](const std::shared_ptr<CoordinateSystem> &cs, const std::vector<double> &params) -> double {
245  Pyston::Context context;
246  context["coordinate_system"] = cs;
247  return dependent(context, params);
248  };
249 
250  m_parameters[p.first] = std::make_shared<FlexibleModelFittingDependentParameter>(
251  p.first, dependent_func, params);
252  }
253 
254  for (auto& p : getDependency<PythonConfig>().getInterpreter().getConstantModels()) {
255  int value_id = py::extract<int>(p.second.attr("value").attr("id"));
256  m_models[p.first] = std::make_shared<FlexibleModelFittingConstantModel>(m_parameters[value_id]);
257  }
258 
259  for (auto& p : getDependency<PythonConfig>().getInterpreter().getPointSourceModels()) {
260  int x_coord_id = py::extract<int>(p.second.attr("x_coord").attr("id"));
261  int y_coord_id = py::extract<int>(p.second.attr("y_coord").attr("id"));
262  int flux_id = py::extract<int>(p.second.attr("flux").attr("id"));
263  m_models[p.first] = std::make_shared<FlexibleModelFittingPointModel>(
264  m_parameters[x_coord_id], m_parameters[y_coord_id], m_parameters[flux_id]);
265  }
266 
267  for (auto& p : getDependency<PythonConfig>().getInterpreter().getSersicModels()) {
268  int x_coord_id = py::extract<int>(p.second.attr("x_coord").attr("id"));
269  int y_coord_id = py::extract<int>(p.second.attr("y_coord").attr("id"));
270  int flux_id = py::extract<int>(p.second.attr("flux").attr("id"));
271  int effective_radius_id = py::extract<int>(p.second.attr("effective_radius").attr("id"));
272  int aspect_ratio_id = py::extract<int>(p.second.attr("aspect_ratio").attr("id"));
273  int angle_id = py::extract<int>(p.second.attr("angle").attr("id"));
274  int n_id = py::extract<int>(p.second.attr("n").attr("id"));
275  m_models[p.first] = std::make_shared<FlexibleModelFittingSersicModel>(
276  m_parameters[x_coord_id], m_parameters[y_coord_id], m_parameters[flux_id], m_parameters[n_id],
277  m_parameters[effective_radius_id], m_parameters[aspect_ratio_id],
278  m_parameters[angle_id]);
279  }
280 
281  for (auto& p : getDependency<PythonConfig>().getInterpreter().getExponentialModels()) {
282  int x_coord_id = py::extract<int>(p.second.attr("x_coord").attr("id"));
283  int y_coord_id = py::extract<int>(p.second.attr("y_coord").attr("id"));
284  int flux_id = py::extract<int>(p.second.attr("flux").attr("id"));
285  int effective_radius_id = py::extract<int>(p.second.attr("effective_radius").attr("id"));
286  int aspect_ratio_id = py::extract<int>(p.second.attr("aspect_ratio").attr("id"));
287  int angle_id = py::extract<int>(p.second.attr("angle").attr("id"));
288  m_models[p.first] = std::make_shared<FlexibleModelFittingExponentialModel>(
289  m_parameters[x_coord_id], m_parameters[y_coord_id], m_parameters[flux_id],
290  m_parameters[effective_radius_id], m_parameters[aspect_ratio_id], m_parameters[angle_id]);
291  }
292 
293  for (auto& p : getDependency<PythonConfig>().getInterpreter().getDeVaucouleursModels()) {
294  int x_coord_id = py::extract<int>(p.second.attr("x_coord").attr("id"));
295  int y_coord_id = py::extract<int>(p.second.attr("y_coord").attr("id"));
296  int flux_id = py::extract<int>(p.second.attr("flux").attr("id"));
297  int effective_radius_id = py::extract<int>(p.second.attr("effective_radius").attr("id"));
298  int aspect_ratio_id = py::extract<int>(p.second.attr("aspect_ratio").attr("id"));
299  int angle_id = py::extract<int>(p.second.attr("angle").attr("id"));
300  m_models[p.first] = std::make_shared<FlexibleModelFittingDevaucouleursModel>(
301  m_parameters[x_coord_id], m_parameters[y_coord_id], m_parameters[flux_id],
302  m_parameters[effective_radius_id], m_parameters[aspect_ratio_id], m_parameters[angle_id]);
303  }
304 
305 #ifdef WITH_ONNX_MODELS
306  for (auto& p : getDependency<PythonConfig>().getInterpreter().getOnnxModels()) {
307  int x_coord_id = py::extract<int>(p.second.attr("x_coord").attr("id"));
308  int y_coord_id = py::extract<int>(p.second.attr("y_coord").attr("id"));
309  int flux_id = py::extract<int>(p.second.attr("flux").attr("id"));
310  int aspect_ratio_id = py::extract<int>(p.second.attr("aspect_ratio").attr("id"));
311  int angle_id = py::extract<int>(p.second.attr("angle").attr("id"));
312  int scale_id = py::extract<int>(p.second.attr("scale").attr("id"));
313 
315  py::dict parameters = py::extract<py::dict>(p.second.attr("params"));
316  py::list names = parameters.keys();
317  for (int i = 0; i < py::len(names); ++i) {
318  std::string name = py::extract<std::string>(names[i]);
319  params[name] = m_parameters[py::extract<int>(parameters[names[i]].attr("id"))];
320  }
321 
323  py::list models = py::extract<py::list>(p.second.attr("models"));
324  for (int i = 0; i < py::len(models); ++i) {
325  std::string model_filename = py::extract<std::string>(models[i]);
326  onnx_models.emplace_back(std::make_shared<OnnxModel>(model_filename));
327 
328  if (onnx_models.back()->getOutputType() != ONNX_TENSOR_ELEMENT_DATA_TYPE_FLOAT ||
329  onnx_models.back()->getOutputShape().size() != 4 ||
330  onnx_models.back()->getOutputShape()[1] != onnx_models.back()->getOutputShape()[2] ||
331  onnx_models.back()->getOutputShape()[3] != 1)
332  {
333  throw Elements::Exception() << "ONNX models for ModelFitting must output a square array of floats";
334  }
335  }
336 
337  m_models[p.first] = std::make_shared<FlexibleModelFittingOnnxModel>(
338  onnx_models, m_parameters[x_coord_id], m_parameters[y_coord_id], m_parameters[flux_id],
339  m_parameters[aspect_ratio_id], m_parameters[angle_id], m_parameters[scale_id], params);
340  }
341 #else
342  if (getDependency<PythonConfig>().getInterpreter().getOnnxModels().size() > 0) {
343  throw Elements::Exception("Trying to use ONNX models but ONNX support is not available");
344  }
345 #endif
346 
347  for (auto& p : getDependency<PythonConfig>().getInterpreter().getFrameModelsMap()) {
349  for (int x : p.second) {
350  model_list.push_back(m_models[x]);
351  }
352  m_frames.push_back(std::make_shared<FlexibleModelFittingFrame>(p.first, model_list));
353  }
354 
355  for (auto& p : getDependency<PythonConfig>().getInterpreter().getPriors()) {
356  auto& prior = p.second;
357  int param_id = py::extract<int>(prior.attr("param"));
358  auto param = m_parameters[param_id];
359 
360  auto value_func = FunctionFromPython<double(const SourceInterface&)>::get(
361  "Prior mean", expr_builder, prior.attr("value")
362  );
363  auto sigma_func = FunctionFromPython<double(const SourceInterface&)>::get(
364  "Prior sigma", expr_builder, prior.attr("sigma")
365  );
366 
367  m_priors[p.first] = std::make_shared<FlexibleModelFittingPrior>(param, value_func, sigma_func);
368  }
369 
370  m_outputs = getDependency<PythonConfig>().getInterpreter().getModelFittingOutputColumns();
371 
372  auto parameters = getDependency<PythonConfig>().getInterpreter().getModelFittingParams();
373  m_least_squares_engine = py::extract<std::string>(parameters["engine"]);
376  }
377  m_max_iterations = py::extract<int>(parameters["max_iterations"]);
378  m_modified_chi_squared_scale = py::extract<double>(parameters["modified_chi_squared_scale"]);
379  m_use_iterative_fitting = py::extract<bool>(parameters["use_iterative_fitting"]);
380  m_meta_iterations = py::extract<int>(parameters["meta_iterations"]);
381  m_deblend_factor = py::extract<double>(parameters["deblend_factor"]);
382  m_meta_iteration_stop = py::extract<double>(parameters["meta_iteration_stop"]);
383 }
384 
386  return m_parameters;
387 }
388 
390  return m_models;
391 }
392 
394  return m_frames;
395 }
396 
398  return m_priors;
399 }
400 
402  return m_outputs;
403 }
404 
405 }
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
T at(T... args)
T back(T... args)
static Logging getLogger(const std::string &name="")
void warn(const std::string &logMessage)
void info(const std::string &logMessage)
const Exception & log(log4cpp::Priority::Value level, Elements::Logging &logger) const
void registerFunction(const std::string &repr, std::function< Signature > functor)
ExpressionTree< Signature > build(const boost::python::object &pyfunc, BuildParams &&... build_params) const
std::string str() const
const std::vector< std::pair< std::string, std::vector< int > > > & getOutputs() const
std::map< int, std::shared_ptr< FlexibleModelFittingParameter > > m_parameters
const std::vector< std::shared_ptr< FlexibleModelFittingFrame > > & getFrames() const
const std::map< int, std::shared_ptr< FlexibleModelFittingModel > > & getModels() const
const std::map< int, std::shared_ptr< FlexibleModelFittingPrior > > & getPriors() const
const std::map< int, std::shared_ptr< FlexibleModelFittingParameter > > & getParameters() const
std::map< int, std::shared_ptr< FlexibleModelFittingPrior > > m_priors
std::map< int, std::shared_ptr< FlexibleModelFittingModel > > m_models
std::vector< std::pair< std::string, std::vector< int > > > m_outputs
void initialize(const UserValues &args) override
std::vector< std::shared_ptr< FlexibleModelFittingFrame > > m_frames
The SourceInterface is an abstract "source" that has properties attached to it.
T clear(T... args)
T emplace_back(T... args)
T empty(T... args)
constexpr double e
static Elements::Logging logger
std::map< std::string, Attribute > AttributeSet
static auto logger
Definition: WCS.cpp:44
static double image_to_world_delta(const Pyston::Context &context, double x, double y)
static double image_to_world_alpha(const Pyston::Context &context, double x, double y)
T push_back(T... args)
static std::function< double(const Pyston::Context &, const std::vector< double > &)> get(const char *readable, Pyston::ExpressionTreeBuilder &builder, py::object py_func, size_t nparams)
static std::function< double(const SourceInterface &)> get(const char *readable, Pyston::ExpressionTreeBuilder &builder, py::object py_func)
static std::function< double(double, const SourceInterface &)> get(const char *readable, Pyston::ExpressionTreeBuilder &builder, py::object py_func)