Add GSkyModel flux_error() and eflux_error() methods (#3074)

Jürgen Knödlseder
1 parent 19518ba9
Showing 7 changed files with 595 additions and 13 deletions
ChangeLog
NEWS
README.md
doc/source/admin/release_history/2.0.rst
include/GModelSky.hpp
pyext/GModelSky.i
src/model/GModelSky.cpp
-2020-11-07
+2020-11-08
  
         * Version 2.0.0 released
         ========================
@@ -14,6 +14,7 @@
         Implement vectorised response computation (#3203)
         Implement analytical gradients for the radial Gaussian model in a
           stacked analysis (#3203)
+        Add GSkyModel flux_error() and eflux_error() methods (#3074)
         Implement pre-computation cache in GModelSpatialGauss (#1300)
  
  
 New Features and Important Changes in GammaLib 2.0.0
  
-7 November 2020
+8 November 2020
  
  
 1. Introduction
@@ -33,6 +33,10 @@ The following methods have been added:
 - GModelSky::flux(GSkyRegion&, GEnergy&, GEnergy&)
 - GModelSky::eflux(GEnergy&, GEnergy&)
 - GModelSky::eflux(GSkyRegion&, GEnergy&, GEnergy&)
+- GModelSky::flux_error(GEnergy&, GEnergy&)
+- GModelSky::flux_error(GSkyRegion&, GEnergy&, GEnergy&)
+- GModelSky::eflux_error(GEnergy&, GEnergy&)
+- GModelSky::eflux_error(GSkyRegion&, GEnergy&, GEnergy&)
 - GModelSpatialRadial::is_energy_dependent()
 - GModelSpatialRadial::is_time_dependent()
 - GOptimizerLM::npars()
@@ -89,6 +93,9 @@ Implemented GMatrixSparse::row() and GMatrixSparse::add_to_row() methods.
  
 8. Model module
 ---------------
+Added methods to GSkyModel to compute the photon flux and energy flux as
+well as their errors (#3439, #3074).
+
 Added analytical gradients to the GModelSpatialRadialGauss class (#3203).
  
 Added is_energy_dependent() and is_time_dependent() methods to GModelSpatialRadial
 GammaLib information
 ====================
-* Version: 2.0.0.dev (7 November 2020)
+* Version: 2.0.0.dev (8 November 2020)
  
 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
  
 .. _2.0:
  
-GammaLib 2.0 release (7 November 2020)
+GammaLib 2.0 release (8 November 2020)
 ======================================
  
 Overview
@@ -34,6 +34,8 @@ Improvements
 * [`3203 <https://cta-redmine.irap.omp.eu/issues/3203>`_] -
   Implement vectorised response computation and analytical gradients for
   the radial Gaussian model in a stacked analysis
+* [`3074 <https://cta-redmine.irap.omp.eu/issues/3074>`_] -
+  Add ``GSkyModel`` ``flux_error()`` and ``eflux_error()`` methods
 * [`1300 <https://cta-redmine.irap.omp.eu/issues/1300>`_] -
   Implement pre-computation cache in GModelSpatialGauss
  
@@ -183,6 +183,16 @@ public:
     double              eflux(const GSkyRegion& region,
                               const GEnergy&    emin,
                               const GEnergy&    emax) const;
+    double              flux_error(const GEnergy& emin,
+                                   const GEnergy& emax) const;
+    double              flux_error(const GSkyRegion& region,
+                                   const GEnergy&    emin,
+                                   const GEnergy&    emax) const;
+    double              eflux_error(const GEnergy& emin,
+                                    const GEnergy& emax) const;
+    double              eflux_error(const GSkyRegion& region,
+                                    const GEnergy&    emin,
+                                    const GEnergy&    emax) const;
  
 protected:
     // Protected methods
@@ -95,6 +95,16 @@ public:
     double              eflux(const GSkyRegion& region,
                               const GEnergy&    emin,
                               const GEnergy&    emax) const;
+    double              flux_error(const GEnergy& emin,
+                                   const GEnergy& emax) const;
+    double              flux_error(const GSkyRegion& region,
+                                   const GEnergy&    emin,
+                                   const GEnergy&    emax) const;
+    double              eflux_error(const GEnergy& emin,
+                                    const GEnergy& emax) const;
+    double              eflux_error(const GSkyRegion& region,
+                                    const GEnergy&    emin,
+                                    const GEnergy&    emax) const;
 };
  
  
@@ -1048,15 +1048,15 @@ GPhotons GModelSky::mc(const double&amp; area,
  
  
 /***********************************************************************//**
- * @brief Return sky model flux
+ * @brief Return sky model photon flux
  *
  * @param[in] emin Minimum energy.
  * @param[in] emax Maximum energy.
- * @return Sky model flux (cm\f$^{-2}\f$ s\f$^{-2}\f$)
+ * @return Sky model photon flux (cm\f$^{-2}\f$ s\f$^{-1}\f$)
  *
- * Returns the spatially integrate flux in the sky model.
+ * Returns the spatially integrate photon flux in the sky model.
  *
- * @todo Take spatial flux normalisation into account
+ * @todo Take spatial photon flux normalisation into account
  ***************************************************************************/
 double GModelSky::flux(const GEnergy& emin, const GEnergy& emax) const
 {
@@ -1107,14 +1107,14 @@ double GModelSky::flux(const GEnergy&amp; emin, const GEnergy&amp; emax) const
  
  
 /***********************************************************************//**
- * @brief Return sky model flux within sky region
+ * @brief Return sky model photon flux within sky region
  *
  * @param[in] region Sky region.
  * @param[in] emin Minimum energy.
  * @param[in] emax Maximum energy.
- * @return Sky model flux (cm\f$^{-2}\f$ s\f$^{-2}\f$)
+ * @return Sky model photon flux (cm\f$^{-2}\f$ s\f$^{-1}\f$)
  *
- * Returns the flux of the sky model within a given sky region.
+ * Returns the photon flux of the sky model within a given sky region.
  ***************************************************************************/
 double GModelSky::flux(const GSkyRegion& region,
                        const GEnergy&    emin,
@@ -1176,7 +1176,7 @@ double GModelSky::flux(const GSkyRegion&amp; region,
  *
  * @param[in] emin Minimum energy.
  * @param[in] emax Maximum energy.
- * @return Sky model energy flux (erg cm\f$^{-2}\f$ s\f$^{-2}\f$)
+ * @return Sky model energy flux (erg cm\f$^{-2}\f$ s\f$^{-1}\f$)
  *
  * Returns the spatially integrate energy flux in the sky model.
  ***************************************************************************/
@@ -1234,7 +1234,7 @@ double GModelSky::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) const
  * @param[in] region Sky region.
  * @param[in] emin Minimum energy.
  * @param[in] emax Maximum energy.
- * @return Sky model flux (cm\f$^{-2}\f$ s\f$^{-2}\f$)
+ * @return Sky model energy flux (erg cm\f$^{-2}\f$ s\f$^{-1}\f$)
  *
  * Returns the energy flux of the sky model within a given sky region.
  ***************************************************************************/
@@ -1294,6 +1294,558 @@ double GModelSky::eflux(const GSkyRegion&amp; region,
  
  
 /***********************************************************************//**
+ * @brief Return sky model photon flux error
+ *
+ * @param[in] emin Minimum energy.
+ * @param[in] emax Maximum energy.
+ * @return Sky model photon flux error (cm\f$^{-2}\f$ s\f$^{-1}\f$)
+ *
+ * Returns the spatially integrate photon flux error in the sky model. The
+ * photon flux error is computed using Gaussian error propagation
+ *
+ * \f[
+ *    \Delta F_{ph} = \sqrt{
+ *    \sum_i \left( \frac{\partial F_{ph}}{\partial p_i} \right)^2
+ *    \left( \Delta p_i \right)^2 }
+ * \f]
+ *
+ * where
+ * \f$F_{ph}\f$ is the sky model photon flux and
+ * \f$p_i\f$ are the sky model parameters.
+ ***************************************************************************/
+double GModelSky::flux_error(const GEnergy& emin, const GEnergy& emax) const
+{
+    // Initialise flux error
+    double flux_error = 0.0;
+
+    // Loop over all model parameters
+    for (int ipar = 0; ipar < size(); ++ipar) {
+
+        // Get non-const model pointer
+        GModelPar* par = const_cast<GModelPar*>(&(*this)[ipar]);
+
+        // Consider only free parameters
+        if (par->is_free()) {
+
+            // Save current model parameter
+            GModelPar current = *par;
+
+            // Autoscale model parameter
+            par->autoscale();
+
+            // Get actual parameter value
+            double x = par->factor_value();
+
+            // Set fixed step size for computation of derivative.
+            const double step_size = 0.0002;
+            double       h         = step_size;
+
+            // Re-adjust the step-size h in case that the initial step size is
+            // larger than the allowed parameter range
+            if (par->has_min() && par->has_max()) {
+                double par_h = par->factor_max() - par->factor_min();
+                if (par_h < h) {
+                    h = par_h;
+                }
+            }
+
+            // Initialise derivative
+            double derivative = 0.0;
+
+            // Continue only if step size is positive
+            if (h > 0.0) {
+
+                // If we are too close to the minimum boundary use a right
+                // sided difference ...
+                if (par->has_min() && ((x-par->factor_min()) < h)) {
+
+                    // Compute fs1
+                    par->factor_value(x);
+                    double fs1 = flux(emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x+h);
+                    double fs2 = flux(emin, emax);
+
+                    // Compute derivative
+                    derivative = (fs1 - fs2) / h;
+
+                }
+
+                // ... otherwise if we are too close to the maximum boundary use
+                // a left sided difference ...
+                else if (par->has_max() && ((par->factor_max()-x) < h)) {
+
+                    // Compute fs1
+                    par->factor_value(x);
+                    double fs1 = flux(emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x-h);
+                    double fs2 = flux(emin, emax);
+
+                    // Compute derivative
+                    derivative = (fs1 - fs2) / h;
+
+                }
+
+                // ... otherwise use a symmetric difference
+                else {
+
+                    // Compute fs1
+                    par->factor_value(x+h);
+                    double fs1 = flux(emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x-h);
+                    double fs2 = flux(emin, emax);
+
+                    // Compute derivative
+                    derivative = 0.5 * (fs1 - fs2) / h;
+
+                }
+
+            } // endif: step size was positive
+
+            // Add contribution
+            if (derivative > 0.0) {
+                double factor = derivative * par->factor_error();
+                flux_error   += factor * factor;
+            }
+
+            // Restore current model parameter
+            *par = current;
+
+        } // endif: parameter was free
+
+    } // endfor: looped over all parameters
+
+    // Take square root if required
+    if (flux_error > 0.0) {
+        flux_error = std::sqrt(flux_error);
+    }
+
+    // Return flux error
+    return flux_error;
+}
+
+
+/***********************************************************************//**
+ * @brief Return sky model photon flux error within sky region
+ *
+ * @param[in] region Sky region.
+ * @param[in] emin Minimum energy.
+ * @param[in] emax Maximum energy.
+ * @return Sky model photon flux error (cm\f$^{-2}\f$ s\f$^{-1}\f$)
+ *
+ * Returns the spatially integrate photon flux error in the sky model. The
+ * photon flux error is computed using Gaussian error propagation
+ *
+ * \f[
+ *    \Delta F_{ph} = \sqrt{
+ *    \sum_i \left( \frac{\partial F_{ph}}{\partial p_i} \right)^2
+ *    \left( \Delta p_i \right)^2 }
+ * \f]
+ *
+ * where
+ * \f$F_{ph}\f$ is the sky model photon flux and
+ * \f$p_i\f$ are the sky model parameters.
+ ***************************************************************************/
+double GModelSky::flux_error(const GSkyRegion& region,
+                             const GEnergy&    emin,
+                             const GEnergy&    emax) const
+{
+    // Initialise flux error
+    double flux_error = 0.0;
+
+    // Loop over all model parameters
+    for (int ipar = 0; ipar < size(); ++ipar) {
+
+        // Get non-const model pointer
+        GModelPar* par = const_cast<GModelPar*>(&(*this)[ipar]);
+
+        // Consider only free parameters
+        if (par->is_free()) {
+
+            // Save current model parameter
+            GModelPar current = *par;
+
+            // Autoscale model parameter
+            par->autoscale();
+
+            // Get actual parameter value
+            double x = par->factor_value();
+
+            // Set fixed step size for computation of derivative.
+            const double step_size = 0.0002;
+            double       h         = step_size;
+
+            // Re-adjust the step-size h in case that the initial step size is
+            // larger than the allowed parameter range
+            if (par->has_min() && par->has_max()) {
+                double par_h = par->factor_max() - par->factor_min();
+                if (par_h < h) {
+                    h = par_h;
+                }
+            }
+
+            // Initialise derivative
+            double derivative = 0.0;
+
+            // Continue only if step size is positive
+            if (h > 0.0) {
+
+                // If we are too close to the minimum boundary use a right
+                // sided difference ...
+                if (par->has_min() && ((x-par->factor_min()) < h)) {
+
+                    // Compute fs1
+                    par->factor_value(x);
+                    double fs1 = flux(region, emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x+h);
+                    double fs2 = flux(region, emin, emax);
+
+                    // Compute derivative
+                    derivative = (fs1 - fs2) / h;
+
+                }
+
+                // ... otherwise if we are too close to the maximum boundary use
+                // a left sided difference ...
+                else if (par->has_max() && ((par->factor_max()-x) < h)) {
+
+                    // Compute fs1
+                    par->factor_value(x);
+                    double fs1 = flux(region, emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x-h);
+                    double fs2 = flux(region, emin, emax);
+
+                    // Compute derivative
+                    derivative = (fs1 - fs2) / h;
+
+                }
+
+                // ... otherwise use a symmetric difference
+                else {
+
+                    // Compute fs1
+                    par->factor_value(x+h);
+                    double fs1 = flux(region, emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x-h);
+                    double fs2 = flux(region, emin, emax);
+
+                    // Compute derivative
+                    derivative = 0.5 * (fs1 - fs2) / h;
+
+                }
+
+            } // endif: step size was positive
+
+            // Add contribution
+            if (derivative > 0.0) {
+                double factor = derivative * par->factor_error();
+                flux_error   += factor * factor;
+            }
+
+            // Restore current model parameter
+            *par = current;
+
+        } // endif: parameter was free
+
+    } // endfor: looped over all parameters
+
+    // Take square root if required
+    if (flux_error > 0.0) {
+        flux_error = std::sqrt(flux_error);
+    }
+    // Return flux error
+    return flux_error;
+}
+
+
+/***********************************************************************//**
+ * @brief Return sky model energy flux error
+ *
+ * @param[in] emin Minimum energy.
+ * @param[in] emax Maximum energy.
+ * @return Sky model energy flux error (erg cm\f$^{-2}\f$ s\f$^{-1}\f$)
+ *
+ * Returns the spatially integrate energy flux error in the sky model. The
+ * energy flux error is computed using Gaussian error propagation
+ *
+ * \f[
+ *    \Delta F_E = \sqrt{
+ *    \sum_i \left( \frac{\partial F_E}{\partial p_i} \right)^2
+ *    \left( \Delta p_i \right)^2 }
+ * \f]
+ *
+ * where
+ * \f$F_E\f$ is the sky model energy flux and
+ * \f$p_i\f$ are the sky model parameters.
+ ***************************************************************************/
+double GModelSky::eflux_error(const GEnergy& emin, const GEnergy& emax) const
+{
+    // Initialise flux error
+    double eflux_error = 0.0;
+
+    // Loop over all model parameters
+    for (int ipar = 0; ipar < size(); ++ipar) {
+
+        // Get non-const model pointer
+        GModelPar* par = const_cast<GModelPar*>(&(*this)[ipar]);
+
+        // Consider only free parameters
+        if (par->is_free()) {
+
+            // Save current model parameter
+            GModelPar current = *par;
+
+            // Autoscale model parameter
+            par->autoscale();
+
+            // Get actual parameter value
+            double x = par->factor_value();
+
+            // Set fixed step size for computation of derivative.
+            const double step_size = 0.0002;
+            double       h         = step_size;
+
+            // Re-adjust the step-size h in case that the initial step size is
+            // larger than the allowed parameter range
+            if (par->has_min() && par->has_max()) {
+                double par_h = par->factor_max() - par->factor_min();
+                if (par_h < h) {
+                    h = par_h;
+                }
+            }
+
+            // Initialise derivative
+            double derivative = 0.0;
+
+            // Continue only if step size is positive
+            if (h > 0.0) {
+
+                // If we are too close to the minimum boundary use a right
+                // sided difference ...
+                if (par->has_min() && ((x-par->factor_min()) < h)) {
+
+                    // Compute fs1
+                    par->factor_value(x);
+                    double fs1 = eflux(emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x+h);
+                    double fs2 = eflux(emin, emax);
+
+                    // Compute derivative
+                    derivative = (fs1 - fs2) / h;
+
+                }
+
+                // ... otherwise if we are too close to the maximum boundary use
+                // a left sided difference ...
+                else if (par->has_max() && ((par->factor_max()-x) < h)) {
+
+                    // Compute fs1
+                    par->factor_value(x);
+                    double fs1 = eflux(emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x-h);
+                    double fs2 = eflux(emin, emax);
+
+                    // Compute derivative
+                    derivative = (fs1 - fs2) / h;
+
+                }
+
+                // ... otherwise use a symmetric difference
+                else {
+
+                    // Compute fs1
+                    par->factor_value(x+h);
+                    double fs1 = eflux(emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x-h);
+                    double fs2 = eflux(emin, emax);
+
+                    // Compute derivative
+                    derivative = 0.5 * (fs1 - fs2) / h;
+
+                }
+
+            } // endif: step size was positive
+
+            // Add contribution
+            if (derivative > 0.0) {
+                double factor = derivative * par->factor_error();
+                eflux_error  += factor * factor;
+            }
+
+            // Restore current model parameter
+            *par = current;
+
+        } // endif: parameter was free
+
+    } // endfor: looped over all parameters
+
+    // Take square root if required
+    if (eflux_error > 0.0) {
+        eflux_error = std::sqrt(eflux_error);
+    }
+
+    // Return flux error
+    return eflux_error;
+}
+
+
+/***********************************************************************//**
+ * @brief Return sky model energy flux error within sky region
+ *
+ * @param[in] region Sky region.
+ * @param[in] emin Minimum energy.
+ * @param[in] emax Maximum energy.
+ * @return Sky model energy flux error (erg cm\f$^{-2}\f$ s\f$^{-1}\f$)
+ *
+ * Returns the spatially integrate energy flux error in the sky model. The
+ * energy flux error is computed using Gaussian error propagation
+ *
+ * \f[
+ *    \Delta F_E = \sqrt{
+ *    \sum_i \left( \frac{\partial F_E}{\partial p_i} \right)^2
+ *    \left( \Delta p_i \right)^2 }
+ * \f]
+ *
+ * where
+ * \f$F_E\f$ is the sky model energy flux and
+ * \f$p_i\f$ are the sky model parameters.
+ ***************************************************************************/
+double GModelSky::eflux_error(const GSkyRegion& region,
+                              const GEnergy&    emin,
+                              const GEnergy&    emax) const
+{
+    // Initialise flux error
+    double eflux_error = 0.0;
+
+    // Loop over all model parameters
+    for (int ipar = 0; ipar < size(); ++ipar) {
+
+        // Get non-const model pointer
+        GModelPar* par = const_cast<GModelPar*>(&(*this)[ipar]);
+
+        // Consider only free parameters
+        if (par->is_free()) {
+
+            // Save current model parameter
+            GModelPar current = *par;
+
+            // Autoscale model parameter
+            par->autoscale();
+
+            // Get actual parameter value
+            double x = par->factor_value();
+
+            // Set fixed step size for computation of derivative.
+            const double step_size = 0.0002;
+            double       h         = step_size;
+
+            // Re-adjust the step-size h in case that the initial step size is
+            // larger than the allowed parameter range
+            if (par->has_min() && par->has_max()) {
+                double par_h = par->factor_max() - par->factor_min();
+                if (par_h < h) {
+                    h = par_h;
+                }
+            }
+
+            // Initialise derivative
+            double derivative = 0.0;
+
+            // Continue only if step size is positive
+            if (h > 0.0) {
+
+                // If we are too close to the minimum boundary use a right
+                // sided difference ...
+                if (par->has_min() && ((x-par->factor_min()) < h)) {
+
+                    // Compute fs1
+                    par->factor_value(x);
+                    double fs1 = eflux(region, emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x+h);
+                    double fs2 = eflux(region, emin, emax);
+
+                    // Compute derivative
+                    derivative = (fs1 - fs2) / h;
+
+                }
+
+                // ... otherwise if we are too close to the maximum boundary use
+                // a left sided difference ...
+                else if (par->has_max() && ((par->factor_max()-x) < h)) {
+
+                    // Compute fs1
+                    par->factor_value(x);
+                    double fs1 = eflux(region, emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x-h);
+                    double fs2 = eflux(region, emin, emax);
+
+                    // Compute derivative
+                    derivative = (fs1 - fs2) / h;
+
+                }
+
+                // ... otherwise use a symmetric difference
+                else {
+
+                    // Compute fs1
+                    par->factor_value(x+h);
+                    double fs1 = eflux(region, emin, emax);
+
+                    // Compute fs2
+                    par->factor_value(x-h);
+                    double fs2 = eflux(region, emin, emax);
+
+                    // Compute derivative
+                    derivative = 0.5 * (fs1 - fs2) / h;
+
+                }
+
+            } // endif: step size was positive
+
+            // Add contribution
+            if (derivative > 0.0) {
+                double factor = derivative * par->factor_error();
+                eflux_error  += factor * factor;
+            }
+
+            // Restore current model parameter
+            *par = current;
+
+        } // endif: parameter was free
+
+    } // endfor: looped over all parameters
+
+    // Take square root if required
+    if (eflux_error > 0.0) {
+        eflux_error = std::sqrt(eflux_error);
+    }
+    // Return flux error
+    return eflux_error;
+}
+
+
+/***********************************************************************//**
  * @brief Print model information
  *
  * @param[in] chatter Chattiness.