diff --git a/ChangeLog b/ChangeLog index 52b76e7..e253ef0 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,11 @@ -2017-12-01 +2018-01-09 * Version 1.5.0 released ========================= + Add GPha::counts_spectrum() and GPha::backscal_spectrum() methods (#2290) + Add GNdarray sign(), operator*= and operator/= (#2290) + Add GSkyMap::counts() and GSkyMap::flux() methods and clip() friend operator (#2290) Add GSkyMap::array() method to Python interface (#2285) CTA observations handle "statistic" attribute in XML file (#2238) Use CSTAT and CHI2 instead of POISSON and GAUSSIAN for statistic (#2237) diff --git a/NEWS b/NEWS index 6ffd063..06e5aa2 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.5.0 -1 December 2017 +9 January 2018 1. Introduction @@ -45,6 +45,9 @@ The following methods have been added: - GEbounds::operator!=(GEbounds&, GEbounds&) - GSkyMap::smooth(std::string&, double&) - GSkyMap::array(int&) - only for Python interface +- GSkyMap::counts() +- GSkyMap::flux() +- clip(GSkyMap&, double&) - GNdarray::shape(std::vector&) - GPha::operator+=(GPha&) - GPha::operator-=(GPha&) @@ -54,6 +57,8 @@ The following methods have been added: - GPha::areascal(int&) - GPha::backscal(int&, double&) - GPha::backscal(int&) +- GPha::counts_spectrum() +- GPha::backscal_spectrum() - operator+(GPha&, GPha&) - operator-(GPha&, GPha&) - operator*(GPha&, double&) @@ -135,7 +140,8 @@ conversion of any spectral model into a spectral node model (#2178). 9. Numerics module ------------------ -None +Added GNdarray sign() friend operator and operator*= and operator/= +(#2290). 10. Observation module @@ -166,6 +172,9 @@ None 12. Sky module -------------- +Added GSkyMap::counts() and GSkyMap::flux() methods and clip() frield +operators that allows clipping of sky map values (#2290) + Added the GSkyMap::array() method to the Python interface to returns the sky map pixels as a list of lists (#2285) @@ -196,6 +205,9 @@ None 16. Xspec module ---------------- +Added GPha::counts_spectrum() and GPha::backscal_spectrum() methods to access +the GPha content (#2290). + Added arithmetic operators for GPha class to add, subtract, scale and divide spectra, and methods to set and retrieve the area and background scaling factors. diff --git a/doc/source/admin/release_history/1.5.rst b/doc/source/admin/release_history/1.5.rst index e4effcf..55991b0 100644 --- a/doc/source/admin/release_history/1.5.rst +++ b/doc/source/admin/release_history/1.5.rst @@ -30,6 +30,10 @@ Bug fixes Improvements ------------ +* [`2290 `_] - + Add ``GSkyMap::counts()`` and ``GSkyMap::flux()`` methods and ``clip()`` + friend operator, ``GNdarray`` ``sign()``, ``operator*=`` and ``operator/=``, + and ``GPha::counts_spectrum()`` and ``GPha::backscal_spectrum()`` methods * [`2285 `_] - Add ``GSkyMap::array()`` method to Python interface * [`2238 `_] - diff --git a/include/GNdarray.hpp b/include/GNdarray.hpp index 36d1ae9..e115adf 100644 --- a/include/GNdarray.hpp +++ b/include/GNdarray.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GNdarray.hpp - N-dimensional array class * * ----------------------------------------------------------------------- * - * copyright (C) 2016-2017 by Juergen Knoedlseder * + * copyright (C) 2016-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * diff --git a/include/GPha.hpp b/include/GPha.hpp index 86e6d03..554efef 100644 --- a/include/GPha.hpp +++ b/include/GPha.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GPha.hpp - XSPEC Pulse Height Analyzer class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2017 by Juergen Knoedlseder * + * copyright (C) 2013-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -122,7 +122,7 @@ protected: void copy_members(const GPha& pha); void free_members(void); void alloc(const int& size); - + // Protected members mutable GFilename m_filename; //!< Filename of origin std::vector m_counts; //!< Counts data diff --git a/include/GSkyMap.hpp b/include/GSkyMap.hpp index c88fea2..4ddc248 100644 --- a/include/GSkyMap.hpp +++ b/include/GSkyMap.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GSkyMap.hpp - Sky map class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2017 by Juergen Knoedlseder * + * copyright (C) 2010-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -92,7 +92,7 @@ class GSkyMap : public GBase { friend GSkyMap log10(const GSkyMap& map); friend GSkyMap abs(const GSkyMap& map); friend GSkyMap sign(const GSkyMap& map); - friend GSkyMap clip(const GSkyMap& map,const double& thresh); + friend GSkyMap clip(const GSkyMap& map, const double& thresh); public: // Constructors and destructors @@ -158,8 +158,10 @@ public: int pix2inx(const GSkyPixel& pixel) const; int dir2inx(const GSkyDir& dir) const; GSkyPixel dir2pix(const GSkyDir& dir) const; + GNdarray counts(void) const; double flux(const int& index, const int& map = 0) const; double flux(const GSkyPixel& pixel, const int& map = 0) const; + GNdarray flux(void) const; double solidangle(const int& index) const; double solidangle(const GSkyPixel& pixel) const; bool contains(const GSkyDir& dir) const; @@ -171,8 +173,6 @@ public: const double* pixels(void) const; GSkyMap extract(const int& map, const int& nmaps = 1) const; void stack_maps(void); - GNdarray total_counts(void) const; - GNdarray total_flux(void) const; void load(const GFilename& filename); void save(const GFilename& filename, const bool& clobber = false) const; diff --git a/inst/cta/include/GCTAOnOffObservation.hpp b/inst/cta/include/GCTAOnOffObservation.hpp index c3de75f..83d9e2f 100644 --- a/inst/cta/include/GCTAOnOffObservation.hpp +++ b/inst/cta/include/GCTAOnOffObservation.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAOnOffObservation.hpp - CTA On/Off observation class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2017 by Chia-Chun Lu & Christoph Deil * + * copyright (C) 2013-2018 by Chia-Chun Lu & Christoph Deil * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -140,14 +140,6 @@ protected: double N_bgd(const GModels& models, const int& ibin, GVector* grad) const; // Likelihood methods - double wstat_value(double non, - double noff, - double alpha, - double ngam, - double& nonpred, - double& nbgd, - double& dlogLdsky, - double& d2logLdsky2) const; virtual double likelihood_cstat(const GModels& models, GVector* gradient, GMatrixSparse* curvature, @@ -156,6 +148,14 @@ protected: GVector* gradient, GMatrixSparse* curvature, double* npred) const; + virtual double wstat_value(const double& non, + const double& noff, + const double& alpha, + const double& ngam, + double& nonpred, + double& nbgd, + double& dlogLdsky, + double& d2logLdsky2) const; // Energy dispersion integration kernel class edisp_kern : public GFunction { diff --git a/inst/cta/src/GCTAOnOffObservation.cpp b/inst/cta/src/GCTAOnOffObservation.cpp index 0b11d85..047d4b0 100644 --- a/inst/cta/src/GCTAOnOffObservation.cpp +++ b/inst/cta/src/GCTAOnOffObservation.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAOnOffObservation.cpp - CTA On/Off observation class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2017 by Chia-Chun Lu & Christoph Deil * + * copyright (C) 2013-2018 by Chia-Chun Lu & Christoph Deil * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -70,12 +70,6 @@ const GObservationRegistry g_onoff_obs_cta_registry(&g_onoff_obs_cta_seed); #define G_LIKELIHOOD "GCTAOnOffObservation::likelihood(GModels&, "\ "GOptimizerPars&, GMatrixSparse&, "\ "GVector&, double&, double&)" -#define G_LIKELIHOOD_CSTAT "GCTAOnOffObservation::likelihood_cstat("\ - "GModels&, GOptimizerPars&, GMatrixSparse&, "\ - "GVector&, double&, double&)" -#define G_LIKELIHOOD_WSTAT "GCTAOnOffObservation::likelihood_wstat("\ - "GModels&, GOptimizerPars&, GMatrixSparse&, "\ - "GVector&, double&, double&)" #define G_SET "GCTAOnOffObservation::set(GCTAObservation&, GSkyDir&, "\ "GSkyRegionMap&, GSkyRegionMap&)" #define G_COMPUTE_ARF "GCTAOnOffObservation::compute_arf(GCTAObservation&, "\ @@ -86,6 +80,13 @@ const GObservationRegistry g_onoff_obs_cta_registry(&g_onoff_obs_cta_seed); "GCTAObservation&, GSkyRegionMap&, GSkyRegionMap&)" #define G_COMPUTE_RMF "GCTAOnOffObservation::compute_rmf(GCTAObservation&, "\ "GSkyRegionMap&)" +#define G_MODEL_BACKGROUND "GCTAOnOffObservation::model_background(GModels&)" +#define G_LIKELIHOOD_CSTAT "GCTAOnOffObservation::likelihood_cstat("\ + "GModels&, GOptimizerPars&, GMatrixSparse&, "\ + "GVector&, double&, double&)" +#define G_LIKELIHOOD_WSTAT "GCTAOnOffObservation::likelihood_wstat("\ + "GModels&, GOptimizerPars&, GMatrixSparse&, "\ + "GVector&, double&, double&)" /* __ Constants __________________________________________________________ */ const double minmod = 1.0e-100; //!< Minimum model value @@ -1736,36 +1737,36 @@ double GCTAOnOffObservation::N_bgd(const GModels& models, } /*********************************************************************** - * @brief Compute predicted gamma-ray counts Pha for given model + * @brief Compute predicted gamma-ray counts for given model * * @param[in] models Model container. - * @return model gamma Pha. - * - * Returns spectrum of predicted number of source events in the On regions + * @return model Pha with predicted gamma-ray counts. * + * Returns spectrum of predicted number of source events in the On + * regions. ***********************************************************************/ GPha GCTAOnOffObservation::model_gamma(const GModels& models) const { - // Get number of model parameters in model container - int npars = models.npars(); + // Get number of model parameters in model container + int npars = models.npars(); - // Create dummy gradient vectors to provide required argument to N_gamma - GVector dummy_grad(npars); + // Create dummy gradient vectors to provide required argument to N_gamma + GVector dummy_grad(npars); - // Initialise output Pha - GEbounds ereco = m_on_spec.ebounds(); - GPha gammas = GPha(ereco); + // Initialise output spectrum + GEbounds ereco = m_on_spec.ebounds(); + GPha gammas = GPha(ereco); - // Loop over all energy bins - for (int i = 0; i < m_on_spec.size(); ++i) { - double ngam = N_gamma(models, i, &dummy_grad); - gammas.fill(ereco.emean(i),ngam); - } + // Loop over all energy bins + for (int i = 0; i < m_on_spec.size(); ++i) { + gammas.fill(ereco.emean(i), N_gamma(models, i, &dummy_grad)); + } - // Return model coun pha - return gammas; + // Return predicted gamma-ray counts spectrum + return gammas; } + /*********************************************************************** * @brief Compute predicted background counts Pha for given model * @@ -1781,70 +1782,73 @@ GPha GCTAOnOffObservation::model_gamma(const GModels& models) const ***********************************************************************/ GPha GCTAOnOffObservation::model_background(const GModels& models) const { - // Get number of model parameters in model container - int npars = models.npars(); + // Get number of model parameters in model container + int npars = models.npars(); - // Create dummy gradient vectors to provide required argument to N_gamma - GVector dummy_grad(npars); + // Create dummy gradient vectors to provide required argument to N_gamma + GVector dummy_grad(npars); - // Initialise output Pha - GEbounds ereco = m_on_spec.ebounds(); - GPha bgds = GPha(ereco); + // Initialise output spectrum + GEbounds ereco = m_on_spec.ebounds(); + GPha bgds = GPha(ereco); - // Extract statistic for this observation - std::string statistic = gammalib::toupper(this->statistic()); + // Extract statistic for this observation + std::string statistic = gammalib::toupper(this->statistic()); - // Loop over all energy bins - for (int i = 0; i < m_on_spec.size(); ++i) { + // Loop over all energy bins + for (int i = 0; i < m_on_spec.size(); ++i) { - // Initialise variable to store number of background counts - double nbgd = 0.0; + // Initialise variable to store number of background counts + double nbgd = 0.0; - // Choose background evaluation method based on fit statistics + // Choose background evaluation method based on fit statistics - // CSTAT - if ((statistic == "POISSON") || (statistic == "CSTAT")) { - nbgd = N_bgd(models, i, &dummy_grad); - } - - // ... or Poisson statistic with measured background - else if (statistic == "WSTAT") { - // Get number of On and Off counts - double non = m_on_spec[i]; - double noff = m_off_spec[i]; - - // Get background scaling - double alpha = m_on_spec.backscal(i); - - // Get number of gamma and background events (and corresponding - // spectral model gradients) - double ngam = N_gamma(models, i, &dummy_grad); - - // Initialise variables for likelihood computation - double nonpred = 0.0; - double dlogLdsky = 0.0; - double d2logLdsky2 = 0.0; - - // Perform likelihood profiling and derive likelihood value - // The number of background counts is updated to the profile value - double logL = wstat_value(non,noff,alpha,ngam, - nonpred,nbgd,dlogLdsky,d2logLdsky2); - } - - // ... or unsupported - else { - std::string msg = "Invalid statistic \""+statistic+"\" encountered. " - "Either specify \"POISSON\", \"CSTAT\" or " - "\"WSTAT\"."; - throw GException::invalid_value(G_LIKELIHOOD, msg); - } - - // Fill background Pha - bgds.fill(ereco.emean(i),nbgd); - } - - // Return model coun pha - return bgds; + // CSTAT + if ((statistic == "POISSON") || (statistic == "CSTAT")) { + nbgd = N_bgd(models, i, &dummy_grad); + } + + // ... or Poisson statistic with measured background + else if (statistic == "WSTAT") { + + // Get number of On and Off counts + double non = m_on_spec[i]; + double noff = m_off_spec[i]; + + // Get background scaling + double alpha = m_on_spec.backscal(i); + + // Get number of gamma and background events (and corresponding + // spectral model gradients) + double ngam = N_gamma(models, i, &dummy_grad); + + // Initialise variables for likelihood computation + double nonpred = 0.0; + double dlogLdsky = 0.0; + double d2logLdsky2 = 0.0; + + // Perform likelihood profiling and derive likelihood value + // The number of background counts is updated to the profile value + double logL = wstat_value(non, noff, alpha, ngam, nonpred, nbgd, + dlogLdsky,d2logLdsky2); + + } // endelse: WSTAT + + // ... or unsupported + else { + std::string msg = "Invalid statistic \""+statistic+"\" encountered. " + "Either specify \"POISSON\", \"CSTAT\" or " + "\"WSTAT\"."; + throw GException::invalid_value(G_MODEL_BACKGROUND, msg); + } + + // Fill background spectrum + bgds.fill(ereco.emean(i), nbgd); + + } // endfor: looped over energy bins + + // Return model count spectrum + return bgds; } @@ -1951,9 +1955,9 @@ GPha GCTAOnOffObservation::model_background(const GModels& models) const * \f] ***********************************************************************/ double GCTAOnOffObservation::likelihood_cstat(const GModels& models, - GVector* gradient, - GMatrixSparse* curvature, - double* npred) const + GVector* gradient, + GMatrixSparse* curvature, + double* npred) const { // Timing measurement #if defined(G_LIKELIHOOD_DEBUG) @@ -2174,6 +2178,188 @@ double GCTAOnOffObservation::likelihood_cstat(const GModels& models, return value; } + + +/***********************************************************************//** + * @brief Evaluate log-likelihood function for On/Off analysis in the + * case of Poisson signal with measured Poisson background + * + * @param[in] models Models. + * @param[in,out] gradient Pointer to gradients. + * @param[in,out] curvature Pointer to curvature matrix. + * @param[in,out] npred Pointer to Npred value. + * @return Log-likelihood value. + * + * @exception GException::invalid_value + * There are no model parameters. + * + * Computes the log-likelihood value for the On/Off observation. The method + * loops over the energy bins to update the function value, its derivatives + * and the curvature matrix. + * + * The number of On counts \f$n_{\rm on}\f$ and Off counts \f$n_{\rm off}\f$ + * are taken from the On and Off spectra, the expected number of gamma-ray + * events \f$\mu_s\f$ is computed from the spectral models of the relevant + * components in the model container (spatial and temporal components + * are ignored so far). See the N_gamma() method for details about the + * model computations. + * + * The estimated number of background counts \f$\mu_b\f$ is derived based + * on the measurement in the Off region by analytical minimization of + * the Poisson likelihood, i.e., it is treated as a nuisance parameter. + * See the wstat_value function for the derivation. + ***********************************************************************/ +double GCTAOnOffObservation::likelihood_wstat(const GModels& models, + GVector* gradient, + GMatrixSparse* curvature, + double* npred) const +{ + // Timing measurement + #if defined(G_LIKELIHOOD_DEBUG) + #ifdef _OPENMP + double t_start = omp_get_wtime(); + #else + clock_t t_start = clock(); + #endif + #endif + + // Initialise statistics + #if defined(G_LIKELIHOOD_DEBUG) + int n_bins = m_on_spec.size(); + int n_used = 0; + int n_small_model = 0; + int n_zero_data = 0; + double sum_data = 0.0; + double sum_model = 0.0; + double init_npred = *npred; + #endif + + // Initialise log-likelihood value + double value = 0.0; + + // Get number of model parameters in model container + int npars = models.npars(); + + // Create model gradient vectors for sky parameters + GVector sky_grad(npars); + + // Allocate working array + GVector colvar(npars); + + // Check that there is at least one parameter + if (npars > 0) { + + // Loop over all energy bins + for (int i = 0; i < m_on_spec.size(); ++i) { + + // Reinitialize working arrays + for (int j = 0; j < npars; ++j) { + sky_grad[j] = 0.0; + } + + // Get number of On and Off counts + double non = m_on_spec[i]; + double noff = m_off_spec[i]; + + // Get background scaling + double alpha = m_on_spec.backscal(i); + + // Get number of gamma and background events (and corresponding + // spectral model gradients) + double ngam = N_gamma(models, i, &sky_grad); + + // Initialise variables for likelihood computation + double nonpred = 0.0; + double nbgd = 0.0; + double dlogLdsky = 0.0; + double d2logLdsky2 = 0.0; + + // Perform likelihood profiling and derive likelihood value + double logL = wstat_value(non, noff, alpha, ngam, + nonpred, nbgd, dlogLdsky, d2logLdsky2); + + // Update global log-likelihood and Npred + value += logL; + *npred += nonpred; + + // Update statistics + #if defined(G_LIKELIHOOD_DEBUG) + n_used++; + sum_data += non; + sum_model += nonpred; + #endif + + // Loop over all parameters + for (int j = 0; j < npars; ++j) { + + // If spectral model for sky component is non-zero and + // non-infinite then handle sky component gradients and + // second derivatives including at least a sky component ... + if (sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) { + + // Gradient + (*gradient)[j] += dlogLdsky * sky_grad[j]; + + // Hessian (from first-order derivatives only) + for (int k = 0; k < npars; ++k) { + + // If spectral model for sky component is non-zero and + // non-infinite then we have the curvature element + // of a sky component + if (sky_grad[k] != 0.0 && + !gammalib::is_infinite(sky_grad[k])) { + colvar[k] = sky_grad[j] * sky_grad[k] * d2logLdsky2; + } + + // ...... else if spectral model for background component + // or neither sky nor background + else { + colvar[k] = 0.0; + } + + } // endfor: Hessian computation + + // Update matrix + curvature->add_to_column(j, colvar); + + } // endif: spectral model is non-zero and non-infinite + + } // endfor: looped over all parameters for derivatives computation + + } // endfor: looped over energy bins + + } // endif: number of parameters was positive + + // ... else there are no parameters, so throw an exception + else { + std::string msg ="No model parameter for the computation of the " + "likelihood in observation \""+this->name()+ + "\" (ID "+this->id()+").\n"; + throw GException::invalid_value(G_LIKELIHOOD_WSTAT, msg); + } + + // Optionally dump gradient and curvature matrix + #if defined(G_LIKELIHOOD_DEBUG) + std::cout << *gradient << std::endl; + std::cout << *curvature << std::endl; + #endif + + // Timing measurement + #if defined(G_LIKELIHOOD_DEBUG) + #ifdef _OPENMP + double t_elapse = omp_get_wtime()-t_start; + #else + double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC; + #endif + std::cout << "GCTAOnOffObservation::optimizer::likelihood_wstat: CPU usage = " + << t_elapse << " sec" << std::endl; + #endif + + // Return + return value; +} + + /***********************************************************************//** * @brief Evaluate log-likelihood value in energy bin for On/Off analysis * by profiling over the number of background counts @@ -2379,279 +2565,96 @@ double GCTAOnOffObservation::likelihood_cstat(const GModels& models, * overestimated if for some bins \f$n_{\rm on}=0\f$ and/or * \f$n_{\rm off}=0\f$ (i.e. if the special cases are encountered). ***********************************************************************/ - -double GCTAOnOffObservation::wstat_value(double non, - double noff, - double alpha, - double ngam, - double& nonpred, - double& nbgd, - double& dlogLdsky, - double& d2logLdsky2) const +double GCTAOnOffObservation::wstat_value(const double& non, + const double& noff, + const double& alpha, + const double& ngam, + double& nonpred, + double& nbgd, + double& dlogLdsky, + double& d2logLdsky2) const { - // Initialise log-likelihood value - double logL; - - // Precompute some values - double alphap1 = alpha + 1.0; - double alpharat = alpha / alphap1; - - // Calculate number of background events, profile likelihood value - // and likelihood derivatives - - // Special case noff = 0 - if (noff == 0.0) { - - // Case A: non = 0. In this case nbgd < 0 hence we set nbgd = 0 - if (non == 0.0) { - nbgd = 0; - nonpred = ngam; - logL = ngam; - dlogLdsky = 1.0; - d2logLdsky2 = 0.0; - } - - // Case B: nbgd is positive - else if (ngam < non * alpharat) { - nbgd = non / alphap1 - ngam / alpha; - nonpred = ngam + alpha * nbgd; - logL = -ngam / alpha - non * std::log(alpharat); - dlogLdsky = -1.0/alpha; - d2logLdsky2 = 0.0; - } - - // Case C: nbgd is zero or negative, hence set nbgd = 0 - else { - nbgd = 0; - nonpred = ngam; - logL = ngam + non * (std::log(non) - std::log(ngam) - 1.0); - dlogLdsky = 1.0 - non / ngam; - d2logLdsky2 = non / (ngam * ngam); - } - - } // endif: noff = 0 - - // Special case non = 0 - else if (non == 0.0) { - nbgd = 0; - nonpred = ngam + alpharat * noff; - logL = ngam + noff * std::log(alphap1); - dlogLdsky = 1.0; - d2logLdsky2 = 0.0; - } // endif: non = 0 - - // General case - else { - - // Compute log-likelihood value - double alphat2 = 2.0 * alpha; - double C = alpha * (non + noff) - alphap1 * ngam; - double D = std::sqrt(C*C + 4.0 * alpha * alphap1 * noff * ngam); - nbgd = (C + D) / (alphat2 * alphap1); - nonpred = ngam + alpha * nbgd; - logL = ngam + alphap1 * nbgd - non - noff - - non * (std::log(nonpred) - std::log(non)) - - noff * (std::log(nbgd) - std::log(noff)); - - // Compute derivatives - double f0 = alphat2 * noff - C; - double dCds = -alphap1; - double dDds = (C * dCds + 2.0 * alphap1 * alpha * noff) / D; - double dbds = (f0 / D - 1.0) / alphat2; - double d2bds2 = (-dCds / D - f0 / (D*D) * dDds) / alphat2; - double f1 = alphap1 - alpha*non/nonpred - noff/nbgd; - double f2 = nonpred * nonpred; - double dpds = 1.0 + alpha * dbds; - double f3 = non / f2 * dpds; - dlogLdsky = 1.0 - non / nonpred + dbds * f1; - d2logLdsky2 = f3 + d2bds2 * f1 + - dbds * (alpha * f3 + noff / (nbgd*nbgd) * dbds); - - } // endelse: general case - - // Return - return logL; -} - - - -/***********************************************************************//** - * @brief Evaluate log-likelihood function for On/Off analysis in the - * case of Poisson signal with measured Poisson background - * - * @param[in] models Models. - * @param[in,out] gradient Pointer to gradients. - * @param[in,out] curvature Pointer to curvature matrix. - * @param[in,out] npred Pointer to Npred value. - * @return Log-likelihood value. - * - * @exception GException::invalid_value - * There are no model parameters. - * - * Computes the log-likelihood value for the On/Off observation. The method - * loops over the energy bins to update the function value, its derivatives - * and the curvature matrix. - * - * The number of On counts \f$n_{\rm on}\f$ and Off counts \f$n_{\rm off}\f$ - * are taken from the On and Off spectra, the expected number of gamma-ray - * events \f$\mu_s\f$ is computed from the spectral models of the relevant - * components in the model container (spatial and temporal components - * are ignored so far). See the N_gamma() method for details about the - * model computations. - * - * The estimated number of background counts \f$\mu_b\f$ is derived based - * on the measurement in the Off region by analytical minimization of - * the Poisson likelihood, i.e., it is treated as a nuisance parameter. - * See the wstat_value function for the derivation. - ***********************************************************************/ -double GCTAOnOffObservation::likelihood_wstat(const GModels& models, - GVector* gradient, - GMatrixSparse* curvature, - double* npred) const -{ - // Timing measurement - #if defined(G_LIKELIHOOD_DEBUG) - #ifdef _OPENMP - double t_start = omp_get_wtime(); - #else - clock_t t_start = clock(); - #endif - #endif - - // Initialise statistics - #if defined(G_LIKELIHOOD_DEBUG) - int n_bins = m_on_spec.size(); - int n_used = 0; - int n_small_model = 0; - int n_zero_data = 0; - double sum_data = 0.0; - double sum_model = 0.0; - double init_npred = *npred; - #endif - // Initialise log-likelihood value - double value = 0.0; - - // Get number of model parameters in model container - int npars = models.npars(); - - // Create model gradient vectors for sky parameters - GVector sky_grad(npars); - - // Allocate working array - GVector colvar(npars); - - // Check that there is at least one parameter - if (npars > 0) { - - // Loop over all energy bins - for (int i = 0; i < m_on_spec.size(); ++i) { - - // Reinitialize working arrays - for (int j = 0; j < npars; ++j) { - sky_grad[j] = 0.0; - } - - // Get number of On and Off counts - double non = m_on_spec[i]; - double noff = m_off_spec[i]; - - // Get background scaling - double alpha = m_on_spec.backscal(i); - - // Get number of gamma and background events (and corresponding - // spectral model gradients) - double ngam = N_gamma(models, i, &sky_grad); - - // Initialise variables for likelihood computation - double nonpred = 0.0; - double nbgd = 0.0; - double dlogLdsky = 0.0; - double d2logLdsky2 = 0.0; - - // Perform likelihood profiling and derive likelihood value - double logL = wstat_value(non,noff,alpha,ngam, - nonpred,nbgd,dlogLdsky,d2logLdsky2); - - // Update global log-likelihood and Npred - value += logL; - *npred += nonpred; - - // Update statistics - #if defined(G_LIKELIHOOD_DEBUG) - n_used++; - sum_data += non; - sum_model += nonpred; - #endif - - // Loop over all parameters - for (int j = 0; j < npars; ++j) { - - // If spectral model for sky component is non-zero and - // non-infinite then handle sky component gradients and - // second derivatives including at least a sky component ... - if (sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) { - - // Gradient - (*gradient)[j] += dlogLdsky * sky_grad[j]; - - // Hessian (from first-order derivatives only) - for (int k = 0; k < npars; ++k) { - - // If spectral model for sky component is non-zero and - // non-infinite then we have the curvature element - // of a sky component - if (sky_grad[k] != 0.0 && - !gammalib::is_infinite(sky_grad[k])) { - colvar[k] = sky_grad[j] * sky_grad[k] * d2logLdsky2; - } - - // ...... else if spectral model for background component - // or neither sky nor background - else { - colvar[k] = 0.0; - } - - } // endfor: Hessian computation + double logL; + + // Precompute some values + double alphap1 = alpha + 1.0; + double alpharat = alpha / alphap1; + + // Calculate number of background events, profile likelihood value + // and likelihood derivatives + // Special case noff = 0 + if (noff == 0.0) { + + // Case A: non = 0. In this case nbgd < 0 hence we set nbgd = 0 + if (non == 0.0) { + nbgd = 0; + nonpred = ngam; + logL = ngam; + dlogLdsky = 1.0; + d2logLdsky2 = 0.0; + } - // Update matrix - curvature->add_to_column(j, colvar); + // Case B: nbgd is positive + else if (ngam < non * alpharat) { + nbgd = non / alphap1 - ngam / alpha; + nonpred = ngam + alpha * nbgd; + logL = -ngam / alpha - non * std::log(alpharat); + dlogLdsky = -1.0/alpha; + d2logLdsky2 = 0.0; + } - } // endif: spectral model is non-zero and non-infinite + // Case C: nbgd is zero or negative, hence set nbgd = 0 + else { + nbgd = 0; + nonpred = ngam; + logL = ngam + non * (std::log(non) - std::log(ngam) - 1.0); + dlogLdsky = 1.0 - non / ngam; + d2logLdsky2 = non / (ngam * ngam); + } - } // endfor: looped over all parameters for derivatives computation + } // endif: noff = 0 - } // endfor: looped over energy bins + // Special case non = 0 + else if (non == 0.0) { + nbgd = 0; + nonpred = ngam + alpharat * noff; + logL = ngam + noff * std::log(alphap1); + dlogLdsky = 1.0; + d2logLdsky2 = 0.0; + } // endif: non = 0 - } // endif: number of parameters was positive - - // ... else there are no parameters, so throw an exception + // General case else { - std::string msg ="No model parameter for the computation of the " - "likelihood in observation \""+this->name()+ - "\" (ID "+this->id()+").\n"; - throw GException::invalid_value(G_LIKELIHOOD_WSTAT, msg); - } - - // Optionally dump gradient and curvature matrix - #if defined(G_LIKELIHOOD_DEBUG) - std::cout << *gradient << std::endl; - std::cout << *curvature << std::endl; - #endif - // Timing measurement - #if defined(G_LIKELIHOOD_DEBUG) - #ifdef _OPENMP - double t_elapse = omp_get_wtime()-t_start; - #else - double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC; - #endif - std::cout << "GCTAOnOffObservation::optimizer::likelihood_wstat: CPU usage = " - << t_elapse << " sec" << std::endl; - #endif + // Compute log-likelihood value + double alphat2 = 2.0 * alpha; + double C = alpha * (non + noff) - alphap1 * ngam; + double D = std::sqrt(C*C + 4.0 * alpha * alphap1 * noff * ngam); + nbgd = (C + D) / (alphat2 * alphap1); + nonpred = ngam + alpha * nbgd; + logL = ngam + alphap1 * nbgd - non - noff - + non * (std::log(nonpred) - std::log(non)) - + noff * (std::log(nbgd) - std::log(noff)); + + // Compute derivatives + double f0 = alphat2 * noff - C; + double dCds = -alphap1; + double dDds = (C * dCds + 2.0 * alphap1 * alpha * noff) / D; + double dbds = (f0 / D - 1.0) / alphat2; + double d2bds2 = (-dCds / D - f0 / (D*D) * dDds) / alphat2; + double f1 = alphap1 - alpha*non/nonpred - noff/nbgd; + double f2 = nonpred * nonpred; + double dpds = 1.0 + alpha * dbds; + double f3 = non / f2 * dpds; + dlogLdsky = 1.0 - non / nonpred + dbds * f1; + d2logLdsky2 = f3 + d2bds2 * f1 + + dbds * (alpha * f3 + noff / (nbgd*nbgd) * dbds); + + } // endelse: general case // Return - return value; + return logL; } diff --git a/pyext/GNdarray.i b/pyext/GNdarray.i index 41ae238..bb16643 100755 --- a/pyext/GNdarray.i +++ b/pyext/GNdarray.i @@ -1,7 +1,7 @@ /*************************************************************************** * GNdarray.i - N-dimensional array class * * ----------------------------------------------------------------------- * - * copyright (C) 2016-2017 by Juergen Knoedlseder * + * copyright (C) 2016-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * diff --git a/pyext/GPha.i b/pyext/GPha.i index 971f6c7..c692008 100644 --- a/pyext/GPha.i +++ b/pyext/GPha.i @@ -1,7 +1,7 @@ /*************************************************************************** * GPha.i - XSPEC Pulse Height Analyzer class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2017 by Juergen Knoedlseder * + * copyright (C) 2013-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * diff --git a/pyext/GSkyMap.i b/pyext/GSkyMap.i index 68b5d10..a6578f0 100644 --- a/pyext/GSkyMap.i +++ b/pyext/GSkyMap.i @@ -1,7 +1,7 @@ /*************************************************************************** * GSkyMap.i - Sky map class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2017 by Juergen Knoedlseder * + * copyright (C) 2010-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -164,8 +164,10 @@ public: int pix2inx(const GSkyPixel& pixel) const; int dir2inx(const GSkyDir& dir) const; GSkyPixel dir2pix(const GSkyDir& dir) const; + GNdarray counts(void); double flux(const int& index, const int& map = 0) const; double flux(const GSkyPixel& pixel, const int& map = 0) const; + GNdarray flux(void); double solidangle(const int& index) const; double solidangle(const GSkyPixel& pixel) const; bool contains(const GSkyDir& dir) const; @@ -177,8 +179,6 @@ public: const double* pixels(void) const; GSkyMap extract(const int& map, const int& nmaps = 1) const; void stack_maps(void); - GNdarray total_counts(void); - GNdarray total_flux(void); void load(const GFilename& filename); void save(const GFilename& filename, const bool& clobber = false) const; @@ -259,7 +259,7 @@ public: return log10(*self); } GSkyMap clip(const double& thresh) { - return clip(*self,thresh); + return clip(*self,thresh); } GSkyMap abs() { return abs(*self); diff --git a/src/numerics/GNdarray.cpp b/src/numerics/GNdarray.cpp index 27c0032..6dba602 100644 --- a/src/numerics/GNdarray.cpp +++ b/src/numerics/GNdarray.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GNdarray.cpp - N-dimensional array class * * ----------------------------------------------------------------------- * - * copyright (C) 2016-2017 by Juergen Knoedlseder * + * copyright (C) 2016-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -333,6 +333,7 @@ GNdarray& GNdarray::operator+=(const GNdarray& array) return *this; } + /***********************************************************************//** * @brief Unary multiplication operator * @@ -355,6 +356,7 @@ GNdarray& GNdarray::operator*=(const GNdarray& array) return *this; } + /***********************************************************************//** * @brief Unary division operator * @@ -370,22 +372,24 @@ GNdarray& GNdarray::operator/=(const GNdarray& array) // Divide elements for (int i = 0; i < m_data.size(); ++i) { - // Only proceed if second array is non zero - if (array.m_data[i] != 0.0){ - m_data[i] /= array.m_data[i]; - } - // Otherwise set ratio to 0 - else { - m_data[i] = 0.0; - } - } + + // Only proceed if second array is non zero + if (array.m_data[i] != 0.0){ + m_data[i] /= array.m_data[i]; + } + + // ... otherwise set ratio to 0 + else { + m_data[i] = 0.0; + } + + } // endfor: looped over array elements // Return array return *this; } - /***********************************************************************//** * @brief Unary subtraction operator * @@ -1221,6 +1225,7 @@ GNdarray log(const GNdarray& array) return result; } + /***********************************************************************//** * @brief Computes sign of array elements * @@ -1234,26 +1239,28 @@ GNdarray sign(const GNdarray& array) // Evaluate each array element for (int i = 0; i < array.m_data.size(); ++i) { - - // Get array content - double content = array.m_data[i]; - // Initialise sign - double sign; + // Get array content + double content = array.m_data[i]; + + // Initialise sign + double sign; - // Handle the 3 cases - if (content < 0.0) { - sign = -1.0; - } - else if (content > 0.0) { - sign = +1.0; - } - else { - sign = 0.0; - } + // Handle the 3 cases + if (content < 0.0) { + sign = -1.0; + } + else if (content > 0.0) { + sign = +1.0; + } + else { + sign = 0.0; + } + // Store sign result.m_data[i] = sign; - } + + } // endfor: looped over all array elements // Return array return result; diff --git a/src/sky/GSkyMap.cpp b/src/sky/GSkyMap.cpp index 9ee542d..cc2ee9b 100644 --- a/src/sky/GSkyMap.cpp +++ b/src/sky/GSkyMap.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GSkyMap.cpp - Sky map class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2017 by Juergen Knoedlseder * + * copyright (C) 2010-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -1373,6 +1373,40 @@ GSkyPixel GSkyMap::dir2pix(const GSkyDir& dir) const /***********************************************************************//** + * @brief Returns array with total number of counts for count maps + * + * @return Array of pixel counts sums. + * + * For a set of n count maps, the method returns a 1D array with n entries + * that correspond to the sum of the pixel counts in each map. + ***************************************************************************/ +GNdarray GSkyMap::counts(void) const +{ + // Initialise output array + GNdarray sum_array = GNdarray(m_num_maps); + + // Loop over maps and store sums in output array + for (int i = 0; i < m_num_maps; ++i) { + + // Initialise sum for this map + double sum = 0.0; + + // Compute map sum + for (int k = 0; k < m_num_pixels; ++k) { + sum += (*this)(k,i); + } + + // Store map sum + sum_array(i) = sum; + + } // endfor: looped over all maps + + // Return sum array + return sum_array; +} + + +/***********************************************************************//** * @brief Returns flux in pixel * * @param[in] index Pixel index [0,...,npix()-1]. @@ -1581,6 +1615,41 @@ double GSkyMap::flux(const GSkyPixel& pixel, const int& map) const /***********************************************************************//** + * @brief Returns array with total flux for sky maps + * + * @return Array of pixel flux sums. + * + * For a set of n sky maps, the method returns a 1D array with n entries + * that correspond to the sum of the pixel flux in each map. + ***************************************************************************/ +GNdarray GSkyMap::flux(void) const +{ + // Initialise output array + GNdarray flux_array = GNdarray(m_num_maps); + + // Loop over maps and store sums in output array + for (int i = 0; i < m_num_maps; ++i) { + + // Initialise flux for this map + double flux = 0.0; + + // Compute map flux + for (int k = 0; k < m_num_pixels; ++k) { + double solidangle = this->solidangle(k); + flux += (*this)(k,i) * solidangle; + } + + // Store map flux + flux_array(i) = flux; + + } // endfor: looped over all maps + + // Return flux array + return flux_array; +} + + +/***********************************************************************//** * @brief Returns solid angle of pixel * * @param[in] index Pixel index [0,...,npix()-1]. @@ -1940,58 +2009,6 @@ void GSkyMap::stack_maps(void) return; } -/***********************************************************************//** - * @brief Returns array with total number of counts for count maps - * - * @return Array of pixel counts sums. - * - * For a set of n count maps, the method returns a 1D array with n entries - * that correspond to the sum of the pixel counts in each map. - ***************************************************************************/ -GNdarray GSkyMap::total_counts(void) const -{ - // Initialise output array - GNdarray collapsed_array = GNdarray(m_num_maps); - - // Loop over maps and store sums in output array - for (int i = 0; i < m_num_maps; ++i) { - double sum = 0.0; - for (int k = 0; k < m_num_pixels; ++k) { - sum += (*this)(k,i); - } - collapsed_array(i) = sum; - } - - // Return pixel - return collapsed_array; -} - -/***********************************************************************//** - * @brief Returns array with total flux for sky maps - * - * @return Array of pixel flux sums. - * - * For a set of n sky maps, the method returns a 1D array with n entries - * that correspond to the sum of the pixel flux in each map. - ***************************************************************************/ -GNdarray GSkyMap::total_flux(void) const -{ - // Initialise output array - GNdarray collapsed_array = GNdarray(m_num_maps); - - // Loop over maps and store sums in output array - for (int i = 0; i < m_num_maps; ++i) { - double flux = 0.0; - for (int k = 0; k < m_num_pixels; ++k) { - double solidangle = this->solidangle(k); - flux += (*this)(k,i) * solidangle; - } - collapsed_array(i) = flux; - } - - // Return pixel - return collapsed_array; -} /***********************************************************************//** * @brief Load skymap from FITS file. @@ -3296,44 +3313,6 @@ GSkyMap sqrt(const GSkyMap& map) return result; } -/***********************************************************************//** - * @brief Clips map at given value - * - * @param[in] map Sky map. - * @param[in] thresh Threshold value. - * @return Sky map containing the value in the original map if >= thresh, - * otherwise thresh - ***************************************************************************/ -GSkyMap clip(const GSkyMap& map,const double& thresh) -{ - // Initialise result vector - GSkyMap result(map); - - // Loop over all maps - for (int i = 0; i < map.nmaps(); ++i) { - - // Loop over all bins - for (int j = 0; j < map.npix(); ++j) { - - // Get the content from the bin - double content = map(j,i); - - // Ic content is below threshold, set to thresh - if (content < thresh) { - result(j,i) = thresh; - } - - else { - result(j,i) = content; - } - - } // endfor: looped over all bins - - } // endfor: looped over all maps - - // Return sky map - return result; -} /***********************************************************************//** * @brief Computes the natural logarithm of sky map elements @@ -3489,3 +3468,44 @@ GSkyMap sign(const GSkyMap& map) // Return sky map return result; } + + +/***********************************************************************//** + * @brief Clips map at given value + * + * @param[in] map Sky map. + * @param[in] thresh Threshold value. + * @return Sky map containing the value in the original map if >= thresh, + * otherwise thresh + ***************************************************************************/ +GSkyMap clip(const GSkyMap& map, const double& thresh) +{ + // Initialise result vector + GSkyMap result(map); + + // Loop over all maps + for (int i = 0; i < map.nmaps(); ++i) { + + // Loop over all bins + for (int j = 0; j < map.npix(); ++j) { + + // Get the content from the bin + double content = map(j,i); + + // If content is below threshold, set to thresh + if (content < thresh) { + result(j,i) = thresh; + } + + // ... otherwise keep content + else { + result(j,i) = content; + } + + } // endfor: looped over all bins + + } // endfor: looped over all maps + + // Return sky map + return result; +} diff --git a/src/xspec/GPha.cpp b/src/xspec/GPha.cpp index 299f2a3..b371ea3 100644 --- a/src/xspec/GPha.cpp +++ b/src/xspec/GPha.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GPha.cpp - XSPEC Pulse Height Analyzer class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2017 by Juergen Knoedlseder * + * copyright (C) 2013-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -582,45 +582,51 @@ double GPha::counts(void) const } /***********************************************************************//** - * @brief Number of counts in spectrum as GNdArray + * @brief Get number of counts in spectrum as GNdarray * - * Returns the number of counts in the spectrum as a GNdArray. + * @return GNdarray of number of counts in spectrum. + * + * Returns the number of counts in the spectrum as GNdarray. ***************************************************************************/ GNdarray GPha::counts_spectrum(void) const { // Initialise array - int size = m_counts.size(); + int size = m_counts.size(); GNdarray counts = GNdarray(size); - + // Compute content for (int i = 0; i < size; ++i) { - counts(i) = m_counts[i]; + counts(i) = m_counts[i]; } // Return counts return counts; } + /***********************************************************************//** - * @brief Number of counts in spectrum as GNdArray + * @brief Get background scaling factors as GNdarray + * + * @return GNdarray of background scaling factors. * - * Returns the number of counts in the spectrum as a GNdArray. + * Returns the background scaling factors as GNdarray. ***************************************************************************/ GNdarray GPha::backscal_spectrum(void) const { // Initialise array - int size = m_backscal.size(); + int size = m_backscal.size(); GNdarray alpha = GNdarray(size); - + // Compute content for (int i = 0; i < size; ++i) { - alpha(i) = m_backscal[i]; + alpha(i) = m_backscal[i]; } - // Return counts + // Return background scaling factors return alpha; } + /***********************************************************************//** * @brief Fill spectrum with a value. * diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp index 0de2349..fc20424 100644 --- a/test/test_GSky.cpp +++ b/test/test_GSky.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_GSky.cpp - Test sky module * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2017 by Juergen Knoedlseder * + * copyright (C) 2010-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -1056,17 +1056,17 @@ void TestGSky::test_GSkyMap(void) test_value(map_stacked.nmaps(), 1, "Test stack_maps() method"); // Test total counts computation - GNdarray counts_spectrum = map_src.total_counts(); + GNdarray counts_spectrum = map_src.counts(); double total_counts = 0.0; for (int i = 0; i < counts_spectrum.size(); ++i) { total_counts += counts_spectrum(i); } - test_value(total_counts, total_src, 1.0e-3, "Test total_counts() method"); - test_value(counts_spectrum.size(), map_src.nmaps(), "Test total_counts() method"); + test_value(total_counts, total_src, 1.0e-3, "Test counts() method"); + test_value(counts_spectrum.size(), map_src.nmaps(), "Test counts() method"); // Test total flux computation - GNdarray flux_spectrum = map_src.total_flux(); - test_value(flux_spectrum.size(), map_src.nmaps(), "Test total_flux() method"); + GNdarray flux_spectrum = map_src.flux(); + test_value(flux_spectrum.size(), map_src.nmaps(), "Test flux() method"); // Test map number changing GSkyMap map_more = map_src;