Improved support for fitting of component scales

Domenico Tiziani
1 parent 703a0409
Showing 4 changed files with 38 additions and 13 deletions
include/GModelSpatialComposite.hpp
src/model/GModelSky.cpp
src/model/GModelSpatialComposite.cpp
test/test_GModel.cpp
@@ -81,7 +81,8 @@ public:
     int                  components(void) const;
     void                 append(const GModelSpatial& component,
                                 const std::string&   name = "",
-                                const double&        scale = 1.0);
+                                const double&        scale = 1.0,
+                                const bool&          freeScale = false);
     const GModelSpatial* component(const int& index) const;
     const GModelSpatial* component(const std::string& name) const;
     double               scale(const int& index) const;
@@ -37,6 +37,7 @@
 #include "GModelSpatialDiffuseCube.hpp"
 #include "GModelSpatialRadial.hpp"
 #include "GModelSpatialElliptical.hpp"
+#include "GModelSpatialComposite.hpp"
 #include "GModelSpatialRegistry.hpp"
 #include "GModelSpectralRegistry.hpp"
 #include "GModelTemporalRegistry.hpp"
@@ -1211,6 +1212,11 @@ void GModelSky::set_type(void)
             m_type = "ExtendedSource";
         }
  
+        // ... otherwise, is spatial model a composite source?
+        else if (dynamic_cast<const GModelSpatialComposite*>(m_spatial) != NULL) {
+            m_type = "CompositeSource";
+        }
+
         // ... otherwise we have a diffuse model
         else {
             m_type = "DiffuseSource";
@@ -360,13 +360,21 @@ void GModelSpatialComposite::read(const GXmlElement&amp; xml)
         // Initialise scale
         double scale = 1.0;
  
+        // Initialise freeScale
+        double freeScale = false;
+
         // Get scale value
         if (spec-> has_attribute("scale")) {
             scale = gammalib::todouble(spec->attribute("scale"));
+
+            // Check if scale is fixed
+            if (spec->has_attribute("free_scale")) {
+                freeScale = gammalib::toint(spec->attribute("free_scale"));
+            }
         }
  
         // Append spatial component to container
-        append(*ptr, name, scale);
+        append(*ptr, name, scale, freeScale);
  
     } // endfor: loop over components
  
@@ -416,7 +424,7 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
         // Find XML element with matching name
         GXmlElement *matching_model = NULL;
         for (int k = 0; k < xml.elements("spatialModel"); ++k) {
-            if (xml.element("spatialModel", k)->attribute("name") == m_names[i]) {
+            if (xml.element("spatialModel", k)->attribute("component") == m_names[i]) {
                 matching_model = xml.element("spatialModel", k);
                 break;
             }
@@ -447,8 +455,11 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
             spatial->write(*matching_model);
  
             // Write scale to XML element if needed
-            if (m_scales_fit[i]->value() != 1.0 || matching_model->has_attribute("scale")) {
+            if (m_scales_fit[i]->value() != 1.0 || matching_model->has_attribute("scale")
+                    || m_scales_fit[i]->is_free()) {
             	matching_model->attribute("scale", gammalib::str(m_scales_fit[i]->value()));
+            	matching_model->attribute("free_scale",
+            	        gammalib::str(m_scales_fit[i]->is_free()));
             }
         }
  
@@ -457,9 +468,14 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
             GXmlElement element("spatialModel");
             spatial->write(element);
  
+            // Write component name
+            element.attribute("component", m_names[i]);
+
             // Write scale to XML element if needed
-            if (m_scales_fit[i]->value() != 1.0) {
+            if (m_scales_fit[i]->value() != 1.0 || m_scales_fit[i]->is_free()) {
                 element.attribute("scale", gammalib::str(m_scales_fit[i]->value()));
+                element.attribute("free_scale",
+                                        gammalib::str(m_scales_fit[i]->is_free()));
             }
             xml.append(element);
         }
@@ -482,12 +498,14 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
  * @param[in] component Spatial model component to append.
  * @param[in] name Name of spatial model (can be empty).
  * @param[in] scale Optional spatial model scaling factor.
+ * @param[in] freeScale Free scale parameter (default:false).
  *
  * Appends a spatial component to the composite model
  ***************************************************************************/
 void GModelSpatialComposite::append(const GModelSpatial& component,
                                     const std::string&   name,
-                                    const double&        scale)
+                                    const double&        scale,
+                                    const bool&          freeScale)
 {
     // Append model container
     m_components.push_back(component.clone());
@@ -532,13 +550,13 @@ void GModelSpatialComposite::append(const GModelSpatial&amp; component,
     // Initialise i-th scale parameter
     m_scales_fit.push_back(new GModelPar());
     m_scales_fit[index]->clear();
-    m_scales_fit[index]->name("scales_fit" + gammalib::str(m_components.size()));
+    m_scales_fit[index]->name(component_name + ":scale");
     m_scales_fit[index]->value(scale);
-    m_scales_fit[index]->range(0.001,1000.);
-    if(index==0)
-    	m_scales_fit[index]->fix();
-    else
+    m_scales_fit[index]->range(1e-10,1e13);
+    if(freeScale)
     	m_scales_fit[index]->free();
+    else
+    	m_scales_fit[index]->fix();
     m_scales_fit[index]->scale(1.0);
     m_scales_fit[index]->gradient(0.0);
     m_scales_fit[index]->has_grad(false);  // for the moment use numerical gradient
@@ -1263,13 +1263,13 @@ void TestGModel::test_spatial_composite(void)
     model1.append(GModelSpatialPointSource(dir1));
     model1.append(GModelSpatialRadialGauss(dir1, 0.2), "", 1.5);
     test_value(model1.components(), 2);
-    test_value(model1.size(), 5);
+    test_value(model1.size(), 7);
  
     // Test XML constructor
     GXml         xml(m_xml_model_spatial_composite);
     GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
     GModelSpatialComposite model2(*element);
-    test_value(model2.size(), 5);
+    test_value(model2.size(), 7);
     test_value(model2.type(), "Composite", "Check model type");
     test_value(model2.scale(0), 1.0);
     test_value(model2.scale(1), 3.0);