Justine Devin / gammalib

Commit b4c68064ede777f379740f4097bdf1d17a047806

Authored by Pierrick Martin
1 parent 80d3ab39

Modified/reorganized the code, added solid angle in background events calculation. Still crashes.

Inline Side-by-side
Showing 4 changed files with 121 additions and 100 deletions
inst/cta/include/GCTAOnOffObservation.hpp
  View file @b4c6806
... ... @@ -89,17 +89,15 @@ public:
89 89 double model_on(const GModels& models,
90 90 const GOptimizerPars& pars,
91 91 int ibin,
92   - GVector* mod_grad);
  92 + GVector* mod_grad) const;
93 93 double model_off(const GModels& models,
94 94 const GOptimizerPars& pars,
95 95 int ibin,
96   - GVector* mod_grad);
  96 + GVector* mod_grad) const;
97 97 double likelihood_poisson_onoff(const GModels& models,
98   - const GOptimizerPars& pars,
99 98 GMatrixSparse* curvature,
100 99 GVector* gradient,
101   - double& value,
102   - double& npred);
  100 + double& npred) const;
103 101  
104 102 protected:
105 103 // Protected methods
... ...
inst/cta/pyext/GCTAOnOffObservation.i
  View file @b4c6806
... ... @@ -93,17 +93,15 @@ public:
93 93 double model_on(const GModels& models,
94 94 const GOptimizerPars& pars,
95 95 int ibin,
96   - GVector* mod_grad);
  96 + GVector* mod_grad) const;
97 97 double model_off(const GModels& models,
98 98 const GOptimizerPars& pars,
99 99 int ibin,
100   - GVector* mod_grad);
  100 + GVector* mod_grad) const;
101 101 double likelihood_poisson_onoff(const GModels& models,
102   - const GOptimizerPars& pars,
103 102 GMatrixSparse* curvature,
104 103 GVector* gradient,
105   - double& value,
106   - double& npred);
  104 + double& npred) const;
107 105 };
108 106  
109 107  
... ...
inst/cta/src/GCTAOnOffObservation.cpp
  View file @b4c6806
... ... @@ -669,7 +669,7 @@ std::string GCTAOnOffObservation::print(const GChatter& chatter) const
669 669 double GCTAOnOffObservation::model_on(const GModels& models,
670 670 const GOptimizerPars& pars,
671 671 int ibin,
672   - GVector* mod_grad)
  672 + GVector* mod_grad) const
673 673 {
674 674 // Get number of parameters
675 675 int npars = pars.size();
... ... @@ -706,7 +706,7 @@ double GCTAOnOffObservation::model_on(const GModels& models,
706 706 const GModel* mptr = models[j];
707 707 if (mptr != NULL) {
708 708  
709   - // Continue only if model applies to specific instrument and observation identifier
  709 + // Continue only if model applies to specific instrument and observation identifier
710 710 if (mptr->is_valid(this->instrument(), this->id())) {
711 711  
712 712 // If this model component is a sky component
... ... @@ -726,25 +726,30 @@ double GCTAOnOffObservation::model_on(const GModels& models,
726 726  
727 727 // Spectral component (the useful one)
728 728 speskyptr=skyptr->spectral();
729   -
730   - // Number of gamma events in model
731   - // (Get flux over energy bin in ph/cm2/s and multiply by effective area and time)
732 729 if (speskyptr != NULL) {
  730 +
  731 + // Number of gamma events in model
  732 + // (Get flux over energy bin in ph/cm2/s and multiply by effective area and time)
733 733 ngam += speskyptr->flux(emin,emax)*m_arf[ibin]*m_ontime;
734   - }
735   -
736   - // Gradients
737   - // Evaluate model at current energy (also fills gradients)
738   - double v=speskyptr->eval_gradients(emean,time);
739   - // Loop over spectral model parameters
740   - for (int k = 0; k < speskyptr->size(); ++k) {
741   - const GModelPar sppar=(*speskyptr)[k];
742   - if ((sppar.is_free()) && (i_par < npars)) {
743   - mod_grad[i_par]=sppar.gradient();
744   - i_par++;
745   - }
  734 + // Debug
  735 + std::cout << "Updating number of gamma events..." << std::endl;
746 736  
747   - } // Looped over parameters of the spectral component of the current model
  737 + // Gradients
  738 + // Evaluate model at current energy (also fills gradients)
  739 + double v=speskyptr->eval_gradients(emean,time);
  740 + // Loop over spectral model parameters
  741 + for (int k = 0; k < speskyptr->size(); ++k) {
  742 + GModelPar sppar=(*speskyptr)[k];
  743 + if (sppar.is_free() && i_par < npars) {
  744 + // Debug
  745 + std::cout << "Setting model gradient for parameter " << i_par << "/" << mod_grad->size() << std::endl;
  746 + mod_grad[i_par]=sppar.gradient();
  747 + i_par++;
  748 + }
  749 +
  750 + } // Looped over parameters of the spectral component of the current model
  751 +
  752 + }
748 753  
749 754 // Increase parameter counter for temporal parameter
750 755 temskyptr=skyptr->temporal();
... ... @@ -793,7 +798,7 @@ double GCTAOnOffObservation::model_on(const GModels&amp; models,
793 798 double GCTAOnOffObservation::model_off(const GModels& models,
794 799 const GOptimizerPars& pars,
795 800 int ibin,
796   - GVector* mod_grad)
  801 + GVector* mod_grad) const
797 802 {
798 803 // Get number of parameters
799 804 int npars = pars.size();
... ... @@ -831,7 +836,7 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
831 836 const GModel* mptr = models[j];
832 837 if (mptr != NULL) {
833 838  
834   - // Continue only if model applies to specific instrument and observation identifier
  839 + // Continue only if model applies to specific instrument and observation identifier
835 840 if (mptr->is_valid(this->instrument(), this->id())) {
836 841  
837 842 // If this model component is a background component
... ... @@ -847,26 +852,37 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
847 852 }
848 853  
849 854 // Spectral component (the useful one)
850   - spebgdptr=bgdptr->spectral();
851   -
852   - // Number of gamma events in model
853   - // (Get flux over energy bin in ph/cm2/s and multiply by effective area and time)
854   - if (spebgdptr != NULL) {
855   - nbgd += spebgdptr->flux(emin,emax)*m_offtime;
856   - }
857   -
858   - // Gradients
859   - // Evaluate model at current energy (also fills gradients)
860   - double v=spebgdptr->eval_gradients(emean,time);
861   - // Loop over spectral model parameters
862   - for (int k = 0; k < spebgdptr->size(); ++k) {
863   - const GModelPar sppar=(*spebgdptr)[k];
864   - if ((sppar.is_free()) && (i_par < npars)) {
865   - mod_grad[i_par]=sppar.gradient();
866   - i_par++;
867   - }
  855 + spebgdptr=bgdptr->spectral();
  856 + if (spebgdptr != NULL) {
868 857  
869   - } // Looped over parameters of the spectral component of the current model
  858 + // Compute total solid angle for OFF regions
  859 + double totsolidangle=0.0;
  860 + for (int m = 0; m < m_off_regions.size(); ++m) {
  861 + const GSkyRegion* skyreg=m_off_regions[m];
  862 + totsolidangle += skyreg->solidangle();
  863 + }
  864 + // Number of gamma events in model
  865 + // (Get flux over energy bin in ph/cm2/s and multiply by effective area and time and solid angle)
  866 + nbgd += spebgdptr->flux(emin,emax)*m_offtime*totsolidangle;
  867 + // Debug
  868 + std::cout << "Updating number of background events..." << std::endl;
  869 +
  870 + // Gradients
  871 + // Evaluate model at current energy (also fills gradients)
  872 + double v=spebgdptr->eval_gradients(emean,time);
  873 + // Loop over spectral model parameters
  874 + for (int k = 0; k < spebgdptr->size(); ++k) {
  875 + GModelPar sppar=(*spebgdptr)[k];
  876 + if (sppar.is_free() && i_par < npars) {
  877 + // Debug
  878 + std::cout << "Setting model gradient for parameter " << i_par << "/" << mod_grad->size() << std::endl;
  879 + mod_grad[i_par]=sppar.gradient();
  880 + i_par++;
  881 + }
  882 +
  883 + } // Looped over parameters of the spectral component of the current model
  884 +
  885 + }
870 886  
871 887 // Increase parameter counter for temporal parameter
872 888 tembgdptr=bgdptr->temporal();
... ... @@ -891,25 +907,37 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
891 907  
892 908 // Spectral component (the useful one)
893 909 spebgdptr=bgdptr->spectral();
894   -
895   - // Number of gamma events in model
896   - // (Get flux over energy bin in ph/cm2/s and multiply by effective area and time)
897 910 if (spebgdptr != NULL) {
898   - nbgd += spebgdptr->flux(emin,emax)*m_offtime;
899   - }
900   -
901   - // Gradients
902   - // Evaluate model at current energy (also fills gradients)
903   - double v=spebgdptr->eval_gradients(emean,time);
904   - // Loop over spectral model parameters
905   - for (int k = 0; k < spebgdptr->size(); ++k) {
906   - const GModelPar sppar=(*spebgdptr)[k];
907   - if ((sppar.is_free()) && (i_par < npars)) {
908   - mod_grad[i_par]=sppar.gradient();
909   - i_par++;
  911 +
  912 + // Compute total solid angle for OFF regions
  913 + double totsolidangle=0.0;
  914 + for (int m = 0; m < m_off_regions.size(); ++m) {
  915 + const GSkyRegion* skyreg=m_off_regions[m];
  916 + totsolidangle += skyreg->solidangle();
910 917 }
911 918  
912   - } // Looped over parameters of the spectral component of the current model
  919 + // Number of gamma events in model
  920 + // (Get flux over energy bin in ph/cm2/s and multiply by effective area and time and solid angle)
  921 + nbgd += spebgdptr->flux(emin,emax)*m_offtime*totsolidangle;
  922 + // Debug
  923 + std::cout << "Updating number of background events..." << std::endl;
  924 +
  925 + // Gradients
  926 + // Evaluate model at current energy (also fills gradients)
  927 + double v=spebgdptr->eval_gradients(emean,time);
  928 + // Loop over spectral model parameters
  929 + for (int k = 0; k < spebgdptr->size(); ++k) {
  930 + GModelPar sppar=(*spebgdptr)[k];
  931 + if (sppar.is_free() && i_par < npars) {
  932 + // Debug
  933 + std::cout << "Setting model gradient for parameter " << i_par << "/" << mod_grad->size() << std::endl;
  934 + mod_grad[i_par]=sppar.gradient();
  935 + i_par++;
  936 + }
  937 +
  938 + } // Looped over parameters of the spectral component of the current model
  939 +
  940 + }
913 941  
914 942 // Increase parameter counter for temporal parameter
915 943 tembgdptr=bgdptr->temporal();
... ... @@ -965,11 +993,9 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
965 993 *
966 994 ***********************************************************************/
967 995 double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels& models,
968   - const GOptimizerPars& pars,
969 996 GMatrixSparse* curvature,
970 997 GVector* gradient,
971   - double& value,
972   - double& npred)
  998 + double& npred) const
973 999 {
974 1000 // Timing measurement
975 1001 #if G_EVAL_TIMING
... ... @@ -984,18 +1010,21 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
984 1010 int n_zero_data = 0;
985 1011 double sum_data = 0.0;
986 1012 double sum_model = 0.0;
987   - double init_value = value;
988 1013 double init_npred = npred;
989 1014 #endif
990 1015  
991   - // Get number of parameters
  1016 + // Initialise likelihood value
  1017 + double value = 0.0;
  1018 +
  1019 + // Get parameter array and number of parameters
  1020 + // (Need to get rid of const qualifier before)
  1021 + GModels* mod=const_cast<GModels*>(&models);
  1022 + const GOptimizerPars pars=mod->pars();
992 1023 int npars = pars.size();
993 1024  
994 1025 // Create model gradients array pointers (one for sky parameters, one for background parameters)
995   - GVector* sky_grad= NULL;
996   - GVector* bgd_grad= NULL;
997   - sky_grad = new GVector(npars);
998   - bgd_grad = new GVector(npars);
  1026 + GVector sky_grad(npars);
  1027 + GVector bgd_grad(npars);
999 1028 // Working arrays
1000 1029 double* values = new double[npars];
1001 1030 int* inx = new int[npars];
... ... @@ -1011,11 +1040,11 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1011 1040  
1012 1041 // Reinitialize working arrays
1013 1042 for (int j = 0; j < npars; ++j) {
1014   - (*sky_grad)[j] = 0.0;
1015   - (*bgd_grad)[j] = 0.0;
  1043 + sky_grad[j] = 0.0;
  1044 + bgd_grad[j] = 0.0;
1016 1045 }
1017 1046  
1018   - // Get energy bin bounds
  1047 + // Get energy bin bounds
1019 1048 const GEnergy emin=m_on_spec.ebounds().emin(i);
1020 1049 const GEnergy emax=m_on_spec.ebounds().emax(i);
1021 1050 const GEnergy emean=m_on_spec.ebounds().elogmean(i);
... ... @@ -1028,8 +1057,8 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1028 1057 double nonpred=0.0;
1029 1058  
1030 1059 // Get number of gamma and background events (and corresponding spectral model gradients)
1031   - nbgd=model_off(models,pars,i,bgd_grad);
1032   - ngam=model_on(models,pars,i,sky_grad);
  1060 + nbgd=model_off(models,pars,i,&bgd_grad);
  1061 + ngam=model_on(models,pars,i,&sky_grad);
1033 1062  
1034 1063 // For debug (to be removed later)
1035 1064 std::cout << "Obs " << m_name << " Bin " << i << " Non=" << non << " Noff=" << noff << " Ngam=" << ngam << " Nbgd=" << nbgd << std::endl;
... ... @@ -1059,8 +1088,8 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1059 1088 // (just needed in call to update hessian matrix)
1060 1089 int ndev = 0;
1061 1090 for (int j = 0; j < npars; ++j) {
1062   - if (((*sky_grad)[j] != 0.0 && !gammalib::is_infinite((*sky_grad)[j])) ||
1063   - ((*bgd_grad)[j] != 0.0 && !gammalib::is_infinite((*bgd_grad)[j]))) {
  1091 + if ((sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) ||
  1092 + (bgd_grad[j] != 0.0 && !gammalib::is_infinite(bgd_grad[j]))) {
1064 1093 inx[ndev] = j;
1065 1094 ndev++;
1066 1095 }
... ... @@ -1077,21 +1106,21 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1077 1106 for (int j = 0; j < npars; ++j) {
1078 1107  
1079 1108 // If spectral model for sky component is non-zero and non-infinite
1080   - if ((*sky_grad)[j] != 0.0 && !gammalib::is_infinite((*sky_grad)[j])) {
  1109 + if (sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) {
1081 1110  
1082 1111 // Gradient
1083   - (*gradient)[j] += sky_factor*(*sky_grad)[j];
  1112 + (*gradient)[j] += sky_factor*sky_grad[j];
1084 1113  
1085 1114 // Hessian (from first-order derivatives only)
1086 1115 for (int k = 0; k < npars; ++k) {
1087 1116  
1088 1117 // If spectral model for sky component is non-zero and non-infinite
1089   - if ((*sky_grad)[k] != 0.0 && !gammalib::is_infinite((*sky_grad)[k])) {
1090   - values[k]=(*sky_grad)[j]*(*sky_grad)[k]*fb;
  1118 + if (sky_grad[k] != 0.0 && !gammalib::is_infinite(sky_grad[k])) {
  1119 + values[k]=sky_grad[j]*sky_grad[k]*fb;
1091 1120  
1092 1121 // If spectral model for sky component is non-zero and non-infinite
1093   - } else if ((*bgd_grad)[k] != 0.0 && !gammalib::is_infinite((*bgd_grad)[k])) {
1094   - values[k]=(*sky_grad)[j]*(*bgd_grad)[k]*fc;
  1122 + } else if (bgd_grad[k] != 0.0 && !gammalib::is_infinite(bgd_grad[k])) {
  1123 + values[k]=sky_grad[j]*bgd_grad[k]*fc;
1095 1124  
1096 1125 // ...else neither sky nor background
1097 1126 } else {
... ... @@ -1104,21 +1133,21 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1104 1133 }
1105 1134  
1106 1135 // If spectral model for sky component is non-zero and non-infinite
1107   - } else if ((*bgd_grad)[j] != 0.0 && !gammalib::is_infinite((*bgd_grad)[j])) {
  1136 + } else if (bgd_grad[j] != 0.0 && !gammalib::is_infinite(bgd_grad[j])) {
1108 1137  
1109 1138 // Gradient
1110   - (*gradient)[j] += bgd_factor*(*bgd_grad)[j];
  1139 + (*gradient)[j] += bgd_factor*bgd_grad[j];
1111 1140  
1112 1141 // Hessian (from first-order derivatives only)
1113 1142 for (int k = 0; k < npars; ++k) {
1114 1143  
1115 1144 // If spectral model for sky component is non-zero and non-infinite
1116   - if ((*sky_grad)[k] != 0.0 && !gammalib::is_infinite((*sky_grad)[k])) {
1117   - values[k]=(*bgd_grad)[j]*(*sky_grad)[k]*fc;
  1145 + if (sky_grad[k] != 0.0 && !gammalib::is_infinite(sky_grad[k])) {
  1146 + values[k]=bgd_grad[j]*bgd_grad[k]*fc;
1118 1147  
1119 1148 // If spectral model for sky component is non-zero and non-infinite
1120   - } else if ((*bgd_grad)[k] != 0.0 && !gammalib::is_infinite((*bgd_grad)[k])) {
1121   - values[k]=(*bgd_grad)[j]*(*bgd_grad)[k]*fd;
  1149 + } else if (bgd_grad[k] != 0.0 && !gammalib::is_infinite(bgd_grad[k])) {
  1150 + values[k]=bgd_grad[j]*bgd_grad[k]*fd;
1122 1151  
1123 1152 // ...else neither sky nor background
1124 1153 } else {
... ... @@ -1156,8 +1185,7 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1156 1185 std::cout << "Number of bins with zero data: " << n_zero_data << std::endl;
1157 1186 std::cout << "Sum of data (ON): " << sum_data << std::endl;
1158 1187 std::cout << "Sum of model (ON): " << sum_model << std::endl;
1159   - std::cout << "Initial statistics: " << init_value << std::endl;
1160   - std::cout << "Delta statistics: " << value-init_value << std::endl;
  1188 + std::cout << "Statistics: " << value << std::endl;
1161 1189 #endif
1162 1190  
1163 1191 // Optionally dump gradient and curvature matrix
... ... @@ -1173,7 +1201,7 @@ double GCTAOnOffObservation::likelihood_poisson_onoff(const GModels&amp;
1173 1201 #else
1174 1202 double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;
1175 1203 #endif
1176   - std::cout << "GCTAOnOffObservation::optimizer::poisson_onoff: CPU usage = "
  1204 + std::cout << "GCTAOnOffObservation::optimizer::likelihood_poisson_onoff: CPU usage = "
1177 1205 << t_elapse << " sec" << std::endl;
1178 1206 #endif
1179 1207  
... ...
inst/cta/src/GCTAOnOffObservations_likelihood.cpp
  View file @b4c6806
... ... @@ -202,15 +202,12 @@ void GCTAOnOffObservations::likelihood::eval(const GOptimizerPars&amp; pars)
202 202 m_curvature = new GMatrixSparse(npars,npars);
203 203  
204 204 // Loop over observations in container
205   - double v=0.0;
206 205 for (int i = 0; i < m_this->size(); ++i) {
207 206  
208 207 // Compute likelihood
209   - v += m_this->m_obs[i]->likelihood_poisson_onoff(m_this->models(),
210   - pars,
  208 + m_value += m_this->m_obs[i]->likelihood_poisson_onoff(m_this->models(),
211 209 m_curvature,
212 210 m_gradient,
213   - m_value,
214 211 m_npred);
215 212  
216 213 } // endfor: looped over observations
... ...