Peter Deiml / gammalib

Commit abbc9eb4d33ea18f685bbead98e4fac3a35bce16

Authored by Jürgen Knödlseder
2 parents 87366eb7 8669abe0

Merge branch 'pull-request-CompositeSpectral' into integration

Inline Side-by-side
Showing 2 changed files with 61 additions and 118 deletions
include/GModelSpatialComposite.hpp
  View file @abbc9eb
... ... @@ -99,7 +99,7 @@ protected:
99 99 std::string m_type; //!< Model type
100 100 std::vector<GModelSpatial*> m_components; //!< Components
101 101 std::vector<std::string> m_names; //!< Component names
102   - std::vector<GModelPar*> m_scales; //!< Component scales
  102 + std::vector<GModelPar> m_scales; //!< Component scales
103 103 mutable GSkyRegionCircle m_region; //!< Bounding circle
104 104 };
105 105  
... ...
src/model/GModelSpatialComposite.cpp
  View file @abbc9eb
... ... @@ -242,12 +242,13 @@ double GModelSpatialComposite::eval(const GPhoton&amp; photon,
242 242 {
243 243 // Initialise value
244 244 double value = 0.0;
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[i]->value();
  248 + value += m_components[i]->eval(photon, gradients) *
  249 + m_scales[i].value();
249 250 }
250   -
  251 +
251 252 // Normalise sum by sum of scales
252 253 if (sum_of_scales() > 0) {
253 254 value /= sum_of_scales();
... ... @@ -282,7 +283,7 @@ GSkyDir GModelSpatialComposite::mc(const GEnergy&amp; energy,
282 283 double sum = 0.0;
283 284 int index = 0;
284 285 for (int i = 0; i < m_scales.size(); ++i) {
285   - sum += m_scales[i]->value();
  286 + sum += m_scales[i].value();
286 287 if (random <= sum) {
287 288 index = i;
288 289 break;
... ... @@ -320,7 +321,7 @@ bool GModelSpatialComposite::contains(const GSkyDir&amp; dir,
320 321 // Loop over all components
321 322 for (int i = 0; i < m_components.size(); ++i) {
322 323 if (m_components[i]->contains(dir, margin) &&
323   - m_scales[i]->value() != 0.0) {
  324 + m_scales[i].value() != 0.0) {
324 325 containment = true;
325 326 break;
326 327 }
... ... @@ -440,16 +441,6 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
440 441 continue;
441 442 }
442 443  
443   - // Find XML element with matching name
444   - GXmlElement *matching_model = NULL;
445   - for (int k = 0; k < xml.elements("spatialModel"); ++k) {
446   - if (xml.element("spatialModel", k)->attribute("component") ==
447   - m_names[i]) {
448   - matching_model = xml.element("spatialModel", k);
449   - break;
450   - }
451   - } // endfor: loop over all XML elements
452   -
453 444 // Create copy of the spatial model
454 445 GModelSpatial* spatial = component->clone();
455 446  
... ... @@ -468,59 +459,56 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
468 459 par.name(name);
469 460 }
470 461  
471   - } // endfor looped over all parameters
472   -
473   - // If an XML element exists then use it ...
474   - if (matching_model != NULL) {
475   -
476   - // Write component into existing XML element
477   - spatial->write(*matching_model);
478   -
479   - // Write scale to XML element if needed
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()) {
486   - matching_model->attribute("scale_error",
487   - gammalib::str(m_scales[i]->error()));
488   - matching_model->attribute("scale_min",
489   - gammalib::str(m_scales[i]->min()));
490   - matching_model->attribute("scale_max",
491   - gammalib::str(m_scales[i]->max()));
492   - }
493   - matching_model->attribute("free_scale",
494   - gammalib::str(m_scales[i]->is_free()));
  462 + } // endfor: looped over all parameters
  463 +
  464 + // Find XML element with matching name
  465 + GXmlElement *matching_model = NULL;
  466 + for (int k = 0; k < xml.elements("spatialModel"); ++k) {
  467 + if (xml.element("spatialModel", k)->attribute("component") ==
  468 + m_names[i]) {
  469 + matching_model = xml.element("spatialModel", k);
  470 + break;
495 471 }
496   - }
  472 + } // endfor: loop over all XML elements
497 473  
498   - // ... otherwise create new XML element
499   - else {
  474 + // If no XML element with matching name was found then create a new
  475 + // XML element and append it
  476 + if (matching_model == NULL) {
500 477  
501   - // ...
  478 + // Create XML element
502 479 GXmlElement element("spatialModel");
503   - spatial->write(element);
504 480  
505 481 // Write component name
506 482 element.attribute("component", m_names[i]);
507 483  
508   - // Write scale to XML element if needed
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()) {
513   - element.attribute("scale_error",
514   - gammalib::str(m_scales[i]->error()));
515   - element.attribute("scale_min",
516   - gammalib::str(m_scales[i]->min()));
517   - element.attribute("scale_max",
518   - gammalib::str(m_scales[i]->max()));
519   - }
520   - element.attribute("free_scale",
521   - gammalib::str(m_scales[i]->is_free()));
522   - }
  484 + // Append XML element
523 485 xml.append(element);
  486 +
  487 + // Set pointer to XML element
  488 + matching_model = xml.element("spatialModel",
  489 + xml.elements("spatialModel")-1);
  490 +
  491 + } // endif: XML element created
  492 +
  493 + // Write component into XML element
  494 + spatial->write(*matching_model);
  495 +
  496 + // Write scale to XML element if needed
  497 + if (m_scales[i].value() != 1.0 ||
  498 + m_scales[i].is_free() ||
  499 + matching_model->has_attribute("scale")) {
  500 + matching_model->attribute("scale",
  501 + gammalib::str(m_scales[i].value()));
  502 + if (m_scales[i].is_free()) {
  503 + matching_model->attribute("scale_error",
  504 + gammalib::str(m_scales[i].error()));
  505 + matching_model->attribute("scale_min",
  506 + gammalib::str(m_scales[i].min()));
  507 + matching_model->attribute("scale_max",
  508 + gammalib::str(m_scales[i].max()));
  509 + }
  510 + matching_model->attribute("free_scale",
  511 + gammalib::str(m_scales[i].is_free()));
524 512 }
525 513  
526 514 // Delete clone
... ... @@ -588,35 +576,15 @@ void GModelSpatialComposite::append(const GModelSpatial&amp; component,
588 576 } // endfor: loop over model parameters
589 577  
590 578 // 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
  579 + GModelPar scale(par);
  580 + scale.name(component_name + ":scale");
  581 + scale.scale(1.0);
  582 + scale.gradient(0.0);
  583 + scale.has_grad(false); // for the moment use numerical gradient
596 584  
597 585 // Put scale parameters in stack
598 586 m_scales.push_back(scale);
599   - m_pars.push_back(scale);
600   -
601   - /*
602   - // Initialise i-th scale parameter
603   - m_scales_fit.push_back(new GModelPar());
604   - m_scales_fit[index]->clear();
605   - m_scales_fit[index]->name(component_name + ":scale");
606   - m_scales_fit[index]->value(scale);
607   - m_scales_fit[index]->range(minScale,maxScale);
608   - if(freeScale)
609   - m_scales_fit[index]->free();
610   - else
611   - m_scales_fit[index]->fix();
612   - m_scales_fit[index]->scale(1.0);
613   - m_scales_fit[index]->gradient(0.0);
614   - m_scales_fit[index]->has_grad(false);
615   -
616   - // Set parameter pointer(s)
617   - GModelPar* par = m_scales_fit.operator[](index);
618   - m_pars.push_back(par);
619   - */
  587 + m_pars.push_back(&scale);
620 588  
621 589 // Return
622 590 return;
... ... @@ -701,7 +669,7 @@ double GModelSpatialComposite::scale(const int&amp; index) const
701 669 }
702 670  
703 671 // Return spatial component scale
704   - return m_scales[index]->value();
  672 + return m_scales[index].value();
705 673 }
706 674  
707 675  
... ... @@ -719,7 +687,7 @@ double GModelSpatialComposite::sum_of_scales(void) const
719 687 for (int i = 0; i < m_scales.size(); ++i) {
720 688  
721 689 // Add scale to sum
722   - sum += m_scales[i]->value();
  690 + sum += m_scales[i].value();
723 691 }
724 692  
725 693 // Return sum
... ... @@ -781,6 +749,7 @@ void GModelSpatialComposite::init_members(void)
781 749 m_names.clear();
782 750 m_scales.clear();
783 751 m_region.clear();
  752 + m_pars.clear();
784 753  
785 754 // Return
786 755 return;
... ... @@ -797,13 +766,9 @@ void GModelSpatialComposite::copy_members(const GModelSpatialComposite&amp; model)
797 766 // Copy members
798 767 m_type = model.m_type;
799 768 m_names = model.m_names;
  769 + m_scales = model.m_scales;
800 770 m_region = model.m_region;
801 771  
802   - // Initialise components
803   - m_components.clear();
804   - m_scales.clear();
805   - m_pars.clear();
806   -
807 772 // Copy components
808 773 for (int i = 0; i < model.m_components.size(); ++i) {
809 774  
... ... @@ -821,15 +786,9 @@ void GModelSpatialComposite::copy_members(const GModelSpatialComposite&amp; model)
821 786  
822 787 } // endfor: looped over all components
823 788  
824   - // Copy m_scales
825   - for (int i = 0; i < model.m_scales.size(); ++i) {
826   -
827   - // Clone scales
828   - GModelPar* par = model.m_scales[i]->clone();
829   - m_scales.push_back(par);
830   -
831   - // Append parameter reference
832   - m_pars.push_back(par);
  789 + // Push back any scales parameters
  790 + for (int i = 0; i < m_scales.size(); ++i) {
  791 + m_pars.push_back(&m_scales[i]);
833 792 }
834 793  
835 794 // Return
... ... @@ -853,22 +812,6 @@ void GModelSpatialComposite::free_members(void)
853 812  
854 813 }
855 814  
856   - // Free m_scales attributes
857   - for (int i = 0; i < m_scales.size(); ++i) {
858   -
859   - // Delete component i
860   - if (m_scales[i] != NULL) {
861   - delete m_scales[i];
862   - }
863   -
864   - // Signal free pointer
865   - m_scales[i] = NULL;
866   - }
867   -
868   - // Empty vectors
869   - m_components.clear();
870   - m_scales.clear();
871   -
872 815 // Return
873 816 return;
874 817 }
... ...