Justine Devin / gammalib

Commit 258b6602492934d13cc047a9362523cf0d3f75eb

Authored by Pierrick Martin
1 parent 0c618d04

Debugged GCTAOnOffObservation class

Inline Side-by-side
Showing 2 changed files with 71 additions and 63 deletions
inst/cta/include/GCTAOnOffObservation.hpp
  View file @258b660
... ... @@ -37,6 +37,11 @@
37 37 #include "GCTAEventAtom.hpp"
38 38 #include "GCTAObservation.hpp"
39 39 #include "GSkyRegions.hpp"
  40 +#include "GModelSpatial.hpp"
  41 +#include "GModelSpectral.hpp"
  42 +#include "GModelTemporal.hpp"
  43 +#include "GModelSky.hpp"
  44 +#include "GCTAModelBackground.hpp"
40 45  
41 46  
42 47 /***********************************************************************//**
... ... @@ -83,16 +88,15 @@ public:
83 88 std::string print(const GChatter& chatter = NORMAL) const;
84 89 double model_on(const GOptimizerPars& pars,
85 90 int ibin,
86   - GVector& mod_grad);
  91 + GVector* mod_grad);
87 92 double model_off(const GOptimizerPars& pars,
88 93 int ibin,
89   - GVector& mod_grad);
  94 + GVector* mod_grad);
90 95 void poisson_onoff(const GOptimizerPars& pars,
91 96 GMatrixSparse& covar,
92 97 GVector& gradient,
93 98 double& value,
94   - double& npred,
95   - GVector& wrk_grad);
  99 + double& npred);
96 100  
97 101 protected:
98 102 // Protected methods
... ... @@ -273,7 +277,7 @@ const GRmf& GCTAOnOffObservation::rmf(void) const
273 277 * @return alpha parameter.
274 278 ***************************************************************************/
275 279 inline
276   -const GRmf& GCTAOnOffObservation::alpha(void) const
  280 +const double GCTAOnOffObservation::alpha(void) const
277 281 {
278 282 return m_alpha;
279 283 }
... ... @@ -285,7 +289,7 @@ const GRmf& GCTAOnOffObservation::alpha(void) const
285 289 * @return ON exposure time
286 290 ***************************************************************************/
287 291 inline
288   -const GRmf& GCTAOnOffObservation::ontime(void) const
  292 +const double GCTAOnOffObservation::ontime(void) const
289 293 {
290 294 return m_ontime;
291 295 }
... ... @@ -297,7 +301,7 @@ const GRmf& GCTAOnOffObservation::ontime(void) const
297 301 * @return OFF exposure time.
298 302 ***************************************************************************/
299 303 inline
300   -const GRmf& GCTAOnOffObservation::offtime(void) const
  304 +const double GCTAOnOffObservation::offtime(void) const
301 305 {
302 306 return m_offtime;
303 307 }
... ...
inst/cta/src/GCTAOnOffObservation.cpp
  View file @258b660
... ... @@ -43,6 +43,9 @@
43 43 "GVector&,"\
44 44 "double&,"\
45 45 "double&)"
  46 +/* __ Constants __________________________________________________________ */
  47 +const double minmod = 1.0e-100; //!< Minimum model value
  48 +const double minerr = 1.0e-100; //!< Minimum statistical error
46 49  
47 50 /* __ Macros _____________________________________________________________ */
48 51  
... ... @@ -664,14 +667,14 @@ std::string GCTAOnOffObservation::print(const GChatter&amp; chatter) const
664 667 ***********************************************************************/
665 668 double GCTAOnOffObservation::model_on(const GOptimizerPars& pars,
666 669 int ibin,
667   - GVector& mod_grad)
  670 + GVector* mod_grad)
668 671 {
669 672 // Get number of parameters
670   - int npars = pars.npars();
  673 + int npars = pars.size();
671 674  
672 675 // Initialize variables (vector has 0.0 values)
673 676 double ngam=0.0;
674   - int ipar=0;
  677 + int i_par=0;
675 678 if (mod_grad != NULL) delete mod_grad;
676 679 mod_grad = new GVector(npars);
677 680  
... ... @@ -679,8 +682,8 @@ double GCTAOnOffObservation::model_on(const GOptimizerPars&amp; pars,
679 682 if (ibin < m_on_spec.size()) {
680 683  
681 684 // Create pointer to model array
682   - GOptimizerPars* parptr=&pars;
683   - models=dynamic_cast<const GModels*>(parptr);
  685 + const GOptimizerPars* parptr=&pars;
  686 + const GModels* models=dynamic_cast<const GModels*>(parptr);
684 687  
685 688 // Check that parameter array if of type GModels
686 689 if (models != NULL) {
... ... @@ -697,30 +700,31 @@ double GCTAOnOffObservation::model_on(const GOptimizerPars&amp; pars,
697 700 for (int j = 0; j < models->size(); ++j) {
698 701  
699 702 // Pointers to model components
700   - GModel* skyptr=NULL;
701 703 GModelSpatial* spaskyptr=NULL;
702 704 GModelSpectral* speskyptr=NULL;
703 705 GModelTemporal* temskyptr=NULL;
704 706  
705 707 // Get model pointer. Continue only if pointer is valid
706   - const GModel* mptr = *models[j];
  708 + const GModel* mptr = (*models)[j];
707 709 if (mptr != NULL) {
708 710  
709 711 // Continue only if model applies to specific instrument and observation identifier
710   - if (mptr->isvalid(m_this->instrument(), m_this->id())) {
  712 + if (mptr->is_valid(this->instrument(), this->id())) {
711 713  
712 714 // If this model component is a sky component
713 715 if (dynamic_cast<const GModelSky*>(mptr) != NULL) {
714   -
715   -
716   - // Increase parameter counter for spatial parameter
717   - spaskyptr=mptr->spatial();
  716 +
  717 + // Set pointer to sky
  718 + const GModelSky* skyptr=dynamic_cast<const GModelSky*>(mptr);
  719 +
  720 + // Increase parameter counter for spatial parameter
  721 + spaskyptr=skyptr->spatial();
718 722 if (spaskyptr != NULL) {
719 723 i_par += spaskyptr->size();
720 724 }
721 725  
722 726 // Spectral component (the useful one)
723   - speskyptr=mptr->spectral();
  727 + speskyptr=skyptr->spectral();
724 728  
725 729 // Number of gamma events in model
726 730 // (Get flux over energy bin in ph/cm2/s and multiply by effective area and time)
... ... @@ -730,19 +734,19 @@ double GCTAOnOffObservation::model_on(const GOptimizerPars&amp; pars,
730 734  
731 735 // Gradients
732 736 // Evaluate model at current energy (also fills gradients)
733   - double v=spskyptr->eval_gradients(emean,time);
  737 + double v=speskyptr->eval_gradients(emean,time);
734 738 // Loop over spectral model parameters
735 739 for (int k = 0; k < speskyptr->size(); ++k) {
736   - GModelPar* sppar=speskyptr[k];
737   - if (sppar->isfree()) {
738   - mod_grad[i_par]=sppar->gradient();
  740 + const GModelPar sppar=(*speskyptr)[k];
  741 + if (sppar.is_free()) {
  742 + mod_grad[i_par]=sppar.gradient();
739 743 i_par++;
740 744 }
741 745  
742 746 } // Looped over parameters of the spectral component of the current model
743 747  
744 748 // Increase parameter counter for temporal parameter
745   - temskyptr=mptr->temporal();
  749 + temskyptr=skyptr->temporal();
746 750 if (temskyptr != NULL) {
747 751 i_par += temskyptr->size();
748 752 }
... ... @@ -787,14 +791,14 @@ double GCTAOnOffObservation::model_on(const GOptimizerPars&amp; pars,
787 791 ***********************************************************************/
788 792 double GCTAOnOffObservation::model_off(const GOptimizerPars& pars,
789 793 int ibin,
790   - GVector& mod_grad)
  794 + GVector* mod_grad)
791 795 {
792 796 // Get number of parameters
793   - int npars = pars.npars();
  797 + int npars = pars.size();
794 798  
795 799 // Initialize variables (vector has 0.0 values)
796 800 double nbgd=0.0;
797   - int ipar=0;
  801 + int i_par=0;
798 802 if (mod_grad != NULL) delete mod_grad;
799 803 mod_grad = new GVector(npars);
800 804  
... ... @@ -802,8 +806,8 @@ double GCTAOnOffObservation::model_off(const GOptimizerPars&amp; pars,
802 806 if (ibin < m_off_spec.size()) {
803 807  
804 808 // Create pointer to model array
805   - GOptimizerPars* parptr=&pars;
806   - models=dynamic_cast<const GModels*>(parptr);
  809 + const GOptimizerPars* parptr=&pars;
  810 + const GModels* models=dynamic_cast<const GModels*>(parptr);
807 811  
808 812 // Check that parameter array if of type GModels
809 813 if (models != NULL) {
... ... @@ -820,30 +824,31 @@ double GCTAOnOffObservation::model_off(const GOptimizerPars&amp; pars,
820 824 for (int j = 0; j < models->size(); ++j) {
821 825  
822 826 // Pointers to model components
823   - GModel* bgdptr=NULL;
824 827 GModelSpatial* spabgdptr=NULL;
825 828 GModelSpectral* spebgdptr=NULL;
826 829 GModelTemporal* tembgdptr=NULL;
827 830  
828 831 // Get model pointer. Continue only if pointer is valid
829   - const GModel* mptr = *models[j];
  832 + const GModel* mptr = (*models)[j];
830 833 if (mptr != NULL) {
831 834  
832 835 // Continue only if model applies to specific instrument and observation identifier
833   - if (mptr->isvalid(m_this->instrument(), m_this->id())) {
  836 + if (mptr->is_valid(this->instrument(), this->id())) {
834 837  
835   - // If this model component is a sky component
836   - if (dynamic_cast<const GModelSky*>(mptr) != NULL) {
  838 + // If this model component is a background component
  839 + if (dynamic_cast<const GCTAModelBackground*>(mptr) != NULL) {
837 840  
  841 + // Set pointer to background
  842 + const GCTAModelBackground* bgdptr=dynamic_cast<const GCTAModelBackground*>(mptr);
838 843  
839 844 // Increase parameter counter for spatial parameter
840   - spabgdptr=mptr->spatial();
  845 + spabgdptr=bgdptr->spatial();
841 846 if (spabgdptr != NULL) {
842 847 i_par += spabgdptr->size();
843 848 }
844 849  
845 850 // Spectral component (the useful one)
846   - spebgdptr=mptr->spectral();
  851 + spebgdptr=bgdptr->spectral();
847 852  
848 853 // Number of gamma events in model
849 854 // (Get flux over energy bin in ph/cm2/s and multiply by effective area and time)
... ... @@ -853,19 +858,19 @@ double GCTAOnOffObservation::model_off(const GOptimizerPars&amp; pars,
853 858  
854 859 // Gradients
855 860 // Evaluate model at current energy (also fills gradients)
856   - double v=spbgdptr->eval_gradients(emean,time);
  861 + double v=spebgdptr->eval_gradients(emean,time);
857 862 // Loop over spectral model parameters
858 863 for (int k = 0; k < spebgdptr->size(); ++k) {
859   - GModelPar* sppar=spebgdptr[k];
860   - if (sppar->isfree()) {
861   - mod_grad[i_par]=sppar->gradient();
  864 + const GModelPar sppar=(*spebgdptr)[k];
  865 + if (sppar.is_free()) {
  866 + mod_grad[i_par]=sppar.gradient();
862 867 i_par++;
863 868 }
864 869  
865 870 } // Looped over parameters of the spectral component of the current model
866 871  
867 872 // Increase parameter counter for temporal parameter
868   - tembgdptr=mptr->temporal();
  873 + tembgdptr=bgdptr->temporal();
869 874 if (tembgdptr != NULL) {
870 875 i_par += tembgdptr->size();
871 876 }
... ... @@ -904,7 +909,6 @@ double GCTAOnOffObservation::model_off(const GOptimizerPars&amp; pars,
904 909 * @param[in,out] gradient Gradient.
905 910 * @param[in,out] value Likelihood value.
906 911 * @param[in,out] npred Number of predicted events.
907   - * @param[in,out] wrk_grad Gradient working array.
908 912 *
909 913 * This method computes the log(likelihood) for one observation in the
910 914 * case of an ON-OFF analysis. The method loops over energy bins to
... ... @@ -942,11 +946,11 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
942 946 #endif
943 947  
944 948 // Get number of parameters
945   - int npars = pars.npars();
  949 + int npars = pars.size();
946 950  
947 951 // Create model gradients array (one for sky parameters, one for background parameters)
948   - sky_grad = new GVector(npars);
949   - bgd_grad = new GVector(npars);
  952 + GVector sky_grad(npars);
  953 + GVector bgd_grad(npars);
950 954 // Working arrays
951 955 double* values = new double[npars];
952 956 int* inx = new int[npars];
... ... @@ -955,8 +959,8 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
955 959 GTime time;
956 960  
957 961 // Create pointer to model array
958   - GOptimizerPars* parptr=&pars;
959   - models=dynamic_cast<const GModels*>(parptr);
  962 + const GOptimizerPars* parptr=&pars;
  963 + const GModels* models=dynamic_cast<const GModels*>(parptr);
960 964  
961 965 // Check that parameter array if of type GModels
962 966 if (models != NULL) {
... ... @@ -977,12 +981,12 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
977 981 double nonpred=0.0;
978 982  
979 983 // Get number of gamma and background events (and corresponding spectral model gradients)
980   - ngam=model_on(pars,i,sky_grad)
981   - nbgd=model_off(pars,i,bgd_grad)
  984 + ngam=model_on(pars,i,&sky_grad);
  985 + nbgd=model_off(pars,i,&bgd_grad);
982 986  
983 987 // Skip bin if model is too small (avoids -Inf or NaN gradients)
984 988 nonpred= ngam+m_alpha*nbgd;
985   - if (nbgd <= m_minmod) || (nonpred <= m_minmod) {
  989 + if ((nbgd <= minmod) || (nonpred <= minmod)) {
986 990 #if G_OPT_DEBUG
987 991 n_small_model++;
988 992 #endif
... ... @@ -991,7 +995,7 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
991 995  
992 996 // Now we have all predicted gamma and background events for current energy bin
993 997 // Update the log(likelihood) and predicted number of events
994   - value += -non*log(nonpred)+nonpred-noff*log(nbgd)+nbgd
  998 + value += -non*log(nonpred)+nonpred-noff*log(nbgd)+nbgd;
995 999 npred += nonpred;
996 1000  
997 1001 // Update statistics
... ... @@ -1005,8 +1009,8 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
1005 1009 // (just needed in call to update hessian matrix)
1006 1010 int ndev = 0;
1007 1011 for (int j = 0; j < npars; ++j) {
1008   - if ((sky_grad[j] != 0.0 && !gammalib::isinfinite(sky_grad[j])) ||
1009   - (bgd_grad[j] != 0.0 && !gammalib::isinfinite(bgd_grad[j]))) {
  1012 + if ((sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) ||
  1013 + (bgd_grad[j] != 0.0 && !gammalib::is_infinite(bgd_grad[j]))) {
1010 1014 inx[ndev] = j;
1011 1015 ndev++;
1012 1016 }
... ... @@ -1015,7 +1019,7 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
1015 1019 // Fill derivatives
1016 1020 double fa=non/nonpred;
1017 1021 double fb=fa/nonpred;
1018   - double fc=m_alpha*fb
  1022 + double fc=m_alpha*fb;
1019 1023 double fd=fc*m_alpha+noff/nbgd/nbgd;
1020 1024 double sky_factor=1.0-fa;
1021 1025 double bgd_factor=1.0+m_alpha-m_alpha*fa-noff/nbgd;
... ... @@ -1023,7 +1027,7 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
1023 1027 for (int j = 0; j < npars; ++j) {
1024 1028  
1025 1029 // If spectral model for sky component is non-zero and non-infinite
1026   - if (sky_grad[j] != 0.0 && !gammalib::isinfinite(sky_grad[j])) {
  1030 + if (sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) {
1027 1031  
1028 1032 // Gradient
1029 1033 gradient[j] += sky_factor*sky_grad[j];
... ... @@ -1032,11 +1036,11 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
1032 1036 for (int k = 0; k < npars; ++k) {
1033 1037  
1034 1038 // If spectral model for sky component is non-zero and non-infinite
1035   - if (sky_grad[k] != 0.0 && !gammalib::isinfinite(sky_grad[k])) {
  1039 + if (sky_grad[k] != 0.0 && !gammalib::is_infinite(sky_grad[k])) {
1036 1040 values[k]=sky_grad[j]*sky_grad[k]*fb;
1037 1041  
1038 1042 // If spectral model for sky component is non-zero and non-infinite
1039   - } else if (bgd_grad[k] != 0.0 && !gammalib::isinfinite(bgd_grad[k])) {
  1043 + } else if (bgd_grad[k] != 0.0 && !gammalib::is_infinite(bgd_grad[k])) {
1040 1044 values[k]=sky_grad[j]*bgd_grad[k]*fc;
1041 1045  
1042 1046 // ...else neither sky nor background
... ... @@ -1050,7 +1054,7 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
1050 1054 }
1051 1055  
1052 1056 // If spectral model for sky component is non-zero and non-infinite
1053   - } else if (bgd_grad[j] != 0.0 && !gammalib::isinfinite(bgd_grad[j])) {
  1057 + } else if (bgd_grad[j] != 0.0 && !gammalib::is_infinite(bgd_grad[j])) {
1054 1058  
1055 1059 // Gradient
1056 1060 gradient[j] += bgd_factor*bgd_grad[j];
... ... @@ -1059,11 +1063,11 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
1059 1063 for (int k = 0; k < npars; ++k) {
1060 1064  
1061 1065 // If spectral model for sky component is non-zero and non-infinite
1062   - if (sky_grad[k] != 0.0 && !gammalib::isinfinite(sky_grad[k])) {
  1066 + if (sky_grad[k] != 0.0 && !gammalib::is_infinite(sky_grad[k])) {
1063 1067 values[k]=bgd_grad[j]*sky_grad[k]*fc;
1064 1068  
1065 1069 // If spectral model for sky component is non-zero and non-infinite
1066   - } else if (bgd_grad[k] != 0.0 && !gammalib::isinfinite(bgd_grad[k])) {
  1070 + } else if (bgd_grad[k] != 0.0 && !gammalib::is_infinite(bgd_grad[k])) {
1067 1071 values[k]=bgd_grad[j]*bgd_grad[k]*fd;
1068 1072  
1069 1073 // ...else neither sky nor background
... ... @@ -1085,8 +1089,8 @@ void GCTAOnOffObservation::poisson_onoff(const GOptimizerPars&amp; pars,
1085 1089 // ... else parameter array is not of type GModel so raise exception
1086 1090 } else {
1087 1091 std::string msg ="Expected a model container for the computation of the "
1088   - "likelihood in observation "+*this.name()+
1089   - " (ID "+*this.id()+").\n";
  1092 + "likelihood in observation "+this->name()+
  1093 + " (ID "+this->id()+").\n";
1090 1094 throw GException::bad_type(G_POISSON_ONOFF,msg);
1091 1095 }
1092 1096  
... ...