Stefan Eschbach / gammalib

Commit e270a3212bd2f98059301bc968d6ff8b72b8d86d

Authored by Jürgen Knödlseder
1 parent db94cfe7

Integrate fitting of composite spatial model scaling factors (#1992)

Inline Side-by-side
Showing 6 changed files with 115 additions and 85 deletions
ChangeLog
  View file @e270a32
1   -2017-04-05
  1 +2017-04-07
2 2  
3 3 * Version 1.3.0 released
4 4 =========================
5 5  
  6 + Allow fitting of scale parameter in composite models (#1992)
6 7 Copy shared Python libraries instead of creating a symbolic link (#1987)
7 8 Added FITS image conversion operators (#1986)
8 9 If code originates from Git, always compile the Python bindings (#1985)
... ...
1 1 New Features and Important Changes in GammaLib 1.3.0.dev1
2 2  
3   -5 April 2017
  3 +7 April 2017
4 4  
5 5  
6 6 1. Introduction
... ... @@ -88,6 +88,9 @@ intensity of a node is set to a non-positive value (#990).
88 88 The GModelSpectralSmoothBrokenPlaw class that implements a smoothly broken
89 89 power law spectrum has been added (#1948).
90 90  
  91 +The scale of the composite spatial models can now be fitted by adding the
  92 +free_scale="1" attribute to the XML file (#1992).
  93 +
91 94  
92 95 9. Numerics module
93 96 ------------------
... ...
include/GModelSpatialComposite.hpp
  View file @e270a32
1 1 /***************************************************************************
2 2 * GModelSpatialPointSource.hpp - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Domenico Tiziani *
  4 + * copyright (C) 2016-2017 by Domenico Tiziani *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -81,10 +81,7 @@ public:
81 81 int components(void) const;
82 82 void append(const GModelSpatial& component,
83 83 const std::string& name = "",
84   - const double& scale = 1.0,
85   - const bool& freeScale = false,
86   - const double& minScale = 1.e-10,
87   - const double& maxScale = 1.e13);
  84 + const GModelPar& par = GModelPar("", 1.0));
88 85 const GModelSpatial* component(const int& index) const;
89 86 const GModelSpatial* component(const std::string& name) const;
90 87 double scale(const int& index) const;
... ... @@ -102,7 +99,7 @@ protected:
102 99 std::string m_type; //!< Model type
103 100 std::vector<GModelSpatial*> m_components; //!< Components
104 101 std::vector<std::string> m_names; //!< Component names
105   - std::vector<GModelPar*> m_scales_fit; //!< scales to be fitted
  102 + std::vector<GModelPar*> m_scales; //!< Component scales
106 103 mutable GSkyRegionCircle m_region; //!< Bounding circle
107 104 };
108 105  
... ...
pyext/GModelSpatialComposite.i
  View file @e270a32
1 1 /***************************************************************************
2 2 * GModelSpatialComposite.i - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Domenico Tiziani *
  4 + * copyright (C) 2016-2017 by Domenico Tiziani *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -65,7 +65,7 @@ public:
65 65 int components(void) const;
66 66 void append(const GModelSpatial& component,
67 67 const std::string& name = "",
68   - const double& scale = 1.0);
  68 + const GModelPar& par = GModelPar("", 1.0));
69 69 const GModelSpatial* component(const int& index) const;
70 70 const GModelSpatial* component(const std::string& name) const;
71 71 double scale(const int& index) const;
... ...
src/model/GModelSpatialComposite.cpp
  View file @e270a32
1 1 /***************************************************************************
2 2 * GModelSpatialComposite.cpp - Spatial composite model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Domenico Tiziani *
  4 + * copyright (C) 2016-2017 by Domenico Tiziani *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -47,7 +47,7 @@ const GModelSpatialRegistry g_spatial_comp_registry(&amp;g_spatial_comp_seed);
47 47 #define G_COMPONENT_INDEX "GModelSpatialComposite::component(int&)"
48 48 #define G_COMPONENT_NAME "GModelSpatialComposite::component(std::string&)"
49 49 #define G_APPEND "GModelSpatialComposite::append(GModelSpatial&, "\
50   - "std::string&)"
  50 + "std::string&, GModelPar&)"
51 51  
52 52 /* __ Macros _____________________________________________________________ */
53 53  
... ... @@ -245,7 +245,7 @@ double GModelSpatialComposite::eval(const GPhoton&amp; photon,
245 245  
246 246 // Sum over all components
247 247 for (int i = 0; i < m_components.size(); ++i) {
248   - value += m_components[i]->eval(photon, gradients) * m_scales_fit[i]->value();
  248 + value += m_components[i]->eval(photon, gradients) * m_scales[i]->value();
249 249 }
250 250  
251 251 // Normalise sum by sum of scales
... ... @@ -281,8 +281,8 @@ GSkyDir GModelSpatialComposite::mc(const GEnergy&amp; energy,
281 281 // Find index of chosen model component
282 282 double sum = 0.0;
283 283 int index = 0;
284   - for (int i = 0; i < m_scales_fit.size(); ++i) {
285   - sum += m_scales_fit[i]->value();
  284 + for (int i = 0; i < m_scales.size(); ++i) {
  285 + sum += m_scales[i]->value();
286 286 if (random <= sum) {
287 287 index = i;
288 288 break;
... ... @@ -319,7 +319,8 @@ bool GModelSpatialComposite::contains(const GSkyDir&amp; dir,
319 319  
320 320 // Loop over all components
321 321 for (int i = 0; i < m_components.size(); ++i) {
322   - if (m_components[i]->contains(dir, margin) && m_scales_fit[i]->value() != 0.0) {
  322 + if (m_components[i]->contains(dir, margin) &&
  323 + m_scales[i]->value() != 0.0) {
323 324 containment = true;
324 325 break;
325 326 }
... ... @@ -338,54 +339,61 @@ bool GModelSpatialComposite::contains(const GSkyDir&amp; dir,
338 339 * Reads the spatial information from an XML element.
339 340 ***************************************************************************/
340 341 void GModelSpatialComposite::read(const GXmlElement& xml)
341   -{
  342 +{
342 343 // Get number of spatial components
343 344 int n_comp = xml.elements("spatialModel");
344   -
  345 +
345 346 // Loop over spatial elements
346 347 for (int i = 0; i < n_comp; ++i) {
347   -
  348 +
348 349 // Get spatial XML element
349 350 const GXmlElement* spec = xml.element("spatialModel", i);
350   -
  351 +
351 352 // Initialise a spatial registry object
352 353 GModelSpatialRegistry registry;
353   -
  354 +
354 355 // Read spatial model
355 356 GModelSpatial* ptr = registry.alloc(*spec);
356   -
  357 +
357 358 // Get component attribute from XML file
358 359 std::string name = spec->attribute("component");
359   -
360   - // Initialise scale
361   - double scale = 1.0;
362 360  
363   - // Initialise scale attributes
364   - double freeScale = false;
365   - double minScale = 1.e-10;
366   - double maxScale = 1.e13;
  361 + // Initialise scale parameter
  362 + GModelPar scale("", 1.0);
  363 + scale.range(1.0e-10, 1.0e10);
  364 + scale.fix();
367 365  
368 366 // Get scale value
369   - if (spec-> has_attribute("scale")) {
370   - scale = gammalib::todouble(spec->attribute("scale"));
  367 + if (spec->has_attribute("scale")) {
  368 +
  369 + // Get scale value
  370 + scale.value(gammalib::todouble(spec->attribute("scale")));
371 371  
372   - // Check if scale is fixed
  372 + // If there is a "free_scale" attribute the fix or free the scaling
  373 + // parameter accordingly
373 374 if (spec->has_attribute("free_scale")) {
374   - freeScale = gammalib::toint(spec->attribute("free_scale"));
  375 + if (gammalib::toint(spec->attribute("free_scale")) == 1) {
  376 + scale.free();
  377 + }
  378 + else {
  379 + scale.fix();
  380 + }
375 381 }
376 382  
377   - // Check if scale minimum is specified
  383 + // If there is a "scale_min" attribute set the minimum accordingly
378 384 if (spec->has_attribute("scale_min")) {
379   - minScale = gammalib::todouble(spec->attribute("scale_min"));
  385 + scale.min(gammalib::todouble(spec->attribute("scale_min")));
380 386 }
381   - // Check if scale maximum is specified
  387 +
  388 + // If there is a "scale_max" attribute set the maximum accordingly
382 389 if (spec->has_attribute("scale_max")) {
383   - maxScale = gammalib::todouble(spec->attribute("scale_max"));
  390 + scale.max(gammalib::todouble(spec->attribute("scale_max")));
384 391 }
385   - }
  392 +
  393 + } // endif: there was a scale attribute
386 394  
387 395 // Append spatial component to container
388   - append(*ptr, name, scale, freeScale, minScale, maxScale);
  396 + append(*ptr, name, scale);
389 397  
390 398 } // endfor: loop over components
391 399  
... ... @@ -400,7 +408,7 @@ void GModelSpatialComposite::read(const GXmlElement&amp; xml)
400 408 * @param[in] xml XML element into which model information is written.
401 409 *
402 410 * @exception GException::model_invalid_spatial
403   - * Existing XML element is not of type 'SkyDirFunction'
  411 + * Existing XML element is not of appropriate type
404 412 * @exception GException::model_invalid_parnum
405 413 * Invalid number of model parameters found in XML element.
406 414 * @exception GException::model_invalid_parnames
... ... @@ -435,7 +443,8 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
435 443 // Find XML element with matching name
436 444 GXmlElement *matching_model = NULL;
437 445 for (int k = 0; k < xml.elements("spatialModel"); ++k) {
438   - if (xml.element("spatialModel", k)->attribute("component") == m_names[i]) {
  446 + if (xml.element("spatialModel", k)->attribute("component") ==
  447 + m_names[i]) {
439 448 matching_model = xml.element("spatialModel", k);
440 449 break;
441 450 }
... ... @@ -463,27 +472,33 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
463 472  
464 473 // If an XML element exists then use it ...
465 474 if (matching_model != NULL) {
  475 +
  476 + // Write component into existing XML element
466 477 spatial->write(*matching_model);
467 478  
468 479 // Write scale to XML element if needed
469   - if (m_scales_fit[i]->value() != 1.0 || matching_model->has_attribute("scale")
470   - || m_scales_fit[i]->is_free()) {
471   - matching_model->attribute("scale", gammalib::str(m_scales_fit[i]->value()));
472   - if (m_scales_fit[i]->is_free()) {
  480 + if (m_scales[i]->value() != 1.0 ||
  481 + matching_model->has_attribute("scale") ||
  482 + m_scales[i]->is_free()) {
  483 + matching_model->attribute("scale",
  484 + gammalib::str(m_scales[i]->value()));
  485 + if (m_scales[i]->is_free()) {
473 486 matching_model->attribute("scale_error",
474   - gammalib::str(m_scales_fit[i]->error()));
  487 + gammalib::str(m_scales[i]->error()));
475 488 matching_model->attribute("scale_min",
476   - gammalib::str(m_scales_fit[i]->min()));
  489 + gammalib::str(m_scales[i]->min()));
477 490 matching_model->attribute("scale_max",
478   - gammalib::str(m_scales_fit[i]->max()));
  491 + gammalib::str(m_scales[i]->max()));
479 492 }
480 493 matching_model->attribute("free_scale",
481   - gammalib::str(m_scales_fit[i]->is_free()));
  494 + gammalib::str(m_scales[i]->is_free()));
482 495 }
483 496 }
484 497  
485 498 // ... otherwise create new XML element
486 499 else {
  500 +
  501 + // ...
487 502 GXmlElement element("spatialModel");
488 503 spatial->write(element);
489 504  
... ... @@ -491,18 +506,19 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
491 506 element.attribute("component", m_names[i]);
492 507  
493 508 // Write scale to XML element if needed
494   - if (m_scales_fit[i]->value() != 1.0 || m_scales_fit[i]->is_free()) {
495   - element.attribute("scale", gammalib::str(m_scales_fit[i]->value()));
496   - if (m_scales_fit[i]->is_free()) {
  509 + if (m_scales[i]->value() != 1.0 ||
  510 + m_scales[i]->is_free()) {
  511 + element.attribute("scale", gammalib::str(m_scales[i]->value()));
  512 + if (m_scales[i]->is_free()) {
497 513 element.attribute("scale_error",
498   - gammalib::str(m_scales_fit[i]->error()));
  514 + gammalib::str(m_scales[i]->error()));
499 515 element.attribute("scale_min",
500   - gammalib::str(m_scales_fit[i]->min()));
  516 + gammalib::str(m_scales[i]->min()));
501 517 element.attribute("scale_max",
502   - gammalib::str(m_scales_fit[i]->max()));
  518 + gammalib::str(m_scales[i]->max()));
503 519 }
504 520 element.attribute("free_scale",
505   - gammalib::str(m_scales_fit[i]->is_free()));
  521 + gammalib::str(m_scales[i]->is_free()));
506 522 }
507 523 xml.append(element);
508 524 }
... ... @@ -523,20 +539,14 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
523 539 * @brief Append spatial component
524 540 *
525 541 * @param[in] component Spatial model component to append.
526   - * @param[in] name Name of spatial model (can be empty).
527   - * @param[in] scale Optional spatial model scaling factor.
528   - * @param[in] freeScale Free scale parameter (default:false).
529   - * @param[in] minScale Minimum for scale parameter (default:1.e-10).
530   - * @param[in] maxScale Maximum for scale parameter (default:1.e13).
  542 + * @param[in] name Name of spatial model.
  543 + * @param[in] par Model scaling parameter.
531 544 *
532 545 * Appends a spatial component to the composite model
533 546 ***************************************************************************/
534 547 void GModelSpatialComposite::append(const GModelSpatial& component,
535 548 const std::string& name,
536   - const double& scale,
537   - const bool& freeScale,
538   - const double& minScale,
539   - const double& maxScale)
  549 + const GModelPar& par)
540 550 {
541 551 // Append model container
542 552 m_components.push_back(component.clone());
... ... @@ -577,7 +587,18 @@ void GModelSpatialComposite::append(const GModelSpatial&amp; component,
577 587  
578 588 } // endfor: loop over model parameters
579 589  
  590 + // Set scaling parameter
  591 + GModelPar* scale = new GModelPar(par);
  592 + scale->name(component_name + ":scale");
  593 + scale->scale(1.0);
  594 + scale->gradient(0.0);
  595 + scale->has_grad(false); // for the moment use numerical gradient
  596 +
  597 + // Put scale parameters in stack
  598 + m_scales.push_back(scale);
  599 + m_pars.push_back(scale);
580 600  
  601 + /*
581 602 // Initialise i-th scale parameter
582 603 m_scales_fit.push_back(new GModelPar());
583 604 m_scales_fit[index]->clear();
... ... @@ -590,11 +611,12 @@ void GModelSpatialComposite::append(const GModelSpatial&amp; component,
590 611 m_scales_fit[index]->fix();
591 612 m_scales_fit[index]->scale(1.0);
592 613 m_scales_fit[index]->gradient(0.0);
593   - m_scales_fit[index]->has_grad(false); // for the moment use numerical gradient
  614 + m_scales_fit[index]->has_grad(false);
594 615  
595 616 // Set parameter pointer(s)
596 617 GModelPar* par = m_scales_fit.operator[](index);
597 618 m_pars.push_back(par);
  619 + */
598 620  
599 621 // Return
600 622 return;
... ... @@ -658,6 +680,7 @@ const GModelSpatial* GModelSpatialComposite::component(const std::string&amp; name)
658 680 return m_components[index];
659 681 }
660 682  
  683 +
661 684 /***********************************************************************//**
662 685 * @brief Returns scale of spatial component
663 686 *
... ... @@ -678,7 +701,7 @@ double GModelSpatialComposite::scale(const int&amp; index) const
678 701 }
679 702  
680 703 // Return spatial component scale
681   - return m_scales_fit[index]->value();
  704 + return m_scales[index]->value();
682 705 }
683 706  
684 707  
... ... @@ -693,10 +716,10 @@ double GModelSpatialComposite::sum_of_scales(void) const
693 716 double sum = 0.0;
694 717  
695 718 // Loop over all component scales
696   - for (int i = 0; i < m_scales_fit.size(); ++i) {
  719 + for (int i = 0; i < m_scales.size(); ++i) {
697 720  
698 721 // Add scale to sum
699   - sum += m_scales_fit[i]->value();
  722 + sum += m_scales[i]->value();
700 723 }
701 724  
702 725 // Return sum
... ... @@ -756,7 +779,7 @@ void GModelSpatialComposite::init_members(void)
756 779 // Initialise other members
757 780 m_components.clear();
758 781 m_names.clear();
759   - m_scales_fit.clear();
  782 + m_scales.clear();
760 783 m_region.clear();
761 784  
762 785 // Return
... ... @@ -772,13 +795,13 @@ void GModelSpatialComposite::init_members(void)
772 795 void GModelSpatialComposite::copy_members(const GModelSpatialComposite& model)
773 796 {
774 797 // Copy members
775   - m_type = model.m_type;
776   - m_names = model.m_names;
777   - m_region = model.m_region;
  798 + m_type = model.m_type;
  799 + m_names = model.m_names;
  800 + m_region = model.m_region;
778 801  
779 802 // Initialise components
780 803 m_components.clear();
781   - m_scales_fit.clear();
  804 + m_scales.clear();
782 805 m_pars.clear();
783 806  
784 807 // Copy components
... ... @@ -798,12 +821,12 @@ void GModelSpatialComposite::copy_members(const GModelSpatialComposite&amp; model)
798 821  
799 822 } // endfor: looped over all components
800 823  
801   - // Copy m_scales_fit
802   - for (int i = 0; i < model.m_scales_fit.size(); ++i) {
  824 + // Copy m_scales
  825 + for (int i = 0; i < model.m_scales.size(); ++i) {
803 826  
804 827 // Clone scales
805   - GModelPar* par = model.m_scales_fit[i]->clone();
806   - m_scales_fit.push_back(par);
  828 + GModelPar* par = model.m_scales[i]->clone();
  829 + m_scales.push_back(par);
807 830  
808 831 // Append parameter reference
809 832 m_pars.push_back(par);
... ... @@ -829,17 +852,23 @@ void GModelSpatialComposite::free_members(void)
829 852 m_components[i] = NULL;
830 853  
831 854 }
832   -
833   - // Free m_scales_fit attributes
834   - for (int i = 0; i < m_scales_fit.size(); ++i) {
  855 +
  856 + // Free m_scales attributes
  857 + for (int i = 0; i < m_scales.size(); ++i) {
835 858  
836 859 // Delete component i
837   - if(m_scales_fit[i] != NULL) delete m_scales_fit[i];
  860 + if (m_scales[i] != NULL) {
  861 + delete m_scales[i];
  862 + }
838 863  
839   - //Signal free pointer
840   - m_scales_fit[i] = NULL;
  864 + // Signal free pointer
  865 + m_scales[i] = NULL;
841 866 }
842 867  
  868 + // Empty vectors
  869 + m_components.clear();
  870 + m_scales.clear();
  871 +
843 872 // Return
844 873 return;
845 874 }
... ...
test/test_GModel.cpp
  View file @e270a32
... ... @@ -1261,7 +1261,7 @@ void TestGModel::test_spatial_composite(void)
1261 1261  
1262 1262 // Test append method
1263 1263 model1.append(GModelSpatialPointSource(dir1));
1264   - model1.append(GModelSpatialRadialGauss(dir1, 0.2), "", 1.5);
  1264 + model1.append(GModelSpatialRadialGauss(dir1, 0.2), "", GModelPar("", 1.5));
1265 1265 test_value(model1.components(), 2);
1266 1266 test_value(model1.size(), 7);
1267 1267  
... ...