Implement radius gradient for radial disk model (#3391)

Jürgen Knödlseder
1 parent 97940c94
Showing 4 changed files with 150 additions and 20 deletions
doc/dev/inst/cta/gammalib_cta.tex
include/GModelSpatialRadialDisk.hpp
src/model/GModelSpatialRadialDisk.cpp
test/dev/test_disk_approximation.py
@@ -151,7 +151,7 @@ and
 \end{equation}
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsubsection{Radial disk}
+\paragraph{Radial disk}
  
 The radial disk function $f(\theta)$ is given by
 \begin{equation}
@@ -168,11 +168,16 @@ with the partial derivatives
 \end{equation}
 and
 \begin{equation}
-\frac{\partial f}{\partial \theta_0}  = ? .
+\frac{\partial f}{\partial \theta_0} = \left \{
+   \begin{array}{l l}
+     \displaystyle -2 \pi f^2(\theta) \sin \theta_0 & \mbox{if $\theta \le \theta_0$} \\ \\
+     \displaystyle 0 & \mbox{if $\theta > \theta_0$}
+   \end{array}
+   \right .
 \end{equation}
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsubsection{Radial ring}
+\paragraph{Radial ring}
  
 The radial ring function $f(\theta)$ is given by
 \begin{equation}
@@ -198,7 +203,7 @@ and
 \end{equation}
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsubsection{Radial Gaussian}
+\paragraph{Radial Gaussian}
  
 The radial Gaussian function $f(\theta)$ is given by
 \begin{equation}
@@ -214,7 +219,7 @@ and
 \end{equation}
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsubsection{Radial shell}
+\paragraph{Radial shell}
  
 The shell function on a sphere is given by
 \begin{equation}
@@ -86,16 +86,17 @@ protected:
     void         init_members(void);
     void         copy_members(const GModelSpatialRadialDisk& model);
     void         free_members(void);
-    void         update(void) const;
+    void         update(const bool& gradients) const;
     virtual void set_region(void) const;
  
     // Protected members
     GModelPar m_radius; //!< Disk radius (degrees)
  
     // Cached members used for pre-computations
-    mutable double m_last_radius;   //!< Last disk radius
-    mutable double m_radius_rad;    //!< Radius in radians
-    mutable double m_norm;          //!< Normalization
+    mutable double m_last_radius;     //!< Last disk radius
+    mutable double m_radius_rad;      //!< Radius in radians
+    mutable double m_norm;            //!< Normalization
+    mutable double m_factor_gradient; //!< Gradient factor
 };
  
  
@@ -51,6 +51,7 @@ const GModelSpatialRegistry   g_radial_disk_legacy_registry(&amp;g_radial_disk_legac
 /* __ Macros _____________________________________________________________ */
  
 /* __ Coding definitions _________________________________________________ */
+//#define G_SMOOTH_RADIUS_GRADIENT            //!< Use smooth radius gradient
  
 /* __ Debug definitions __________________________________________________ */
  
@@ -293,11 +294,36 @@ double GModelSpatialRadialDisk::eval(const double&amp;  theta,
                                      const bool&    gradients) const
 {
     // Update precomputation cache
-    update();
+    update(gradients);
  
     // Set value
     double value = (theta <= m_radius_rad) ? m_norm : 0.0;
  
+    // Optionally compute gradients
+    if (gradients) {
+
+        // Alternative radius
+        #if defined(G_SMOOTH_RADIUS_GRADIENT)
+        double g_radius = 0.0;
+        if (m_radius.is_free()) {
+            const double k = 1000.0;
+            double       f1 = std::exp(-k * (m_radius_rad - theta));
+            double       f2 = 1.0 + f1;
+            g_radius        = m_norm * k / (f2 * f2) * m_radius.scale() *
+                              gammalib::deg2rad;
+        }
+
+        // Evaluate radius gradient
+        #else
+        double g_radius = (m_radius.is_free() && theta <= m_radius_rad)
+                          ? m_factor_gradient : 0.0;
+        #endif
+
+        // Set radius gradient
+        m_radius.factor_gradient(g_radius);
+
+    } // endif: gradient computation was requested
+
     // Compile option: Check for NaN/Inf
     #if defined(G_NAN_CHECK)
     if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
@@ -522,16 +548,21 @@ void GModelSpatialRadialDisk::init_members(void)
     m_radius.free();
     m_radius.scale(1.0);
     m_radius.gradient(0.0);
-    m_radius.has_grad(false);  // Radial components never have gradients
+
+    // Signal that this model provides analytical parameter gradients
+    //m_ra.has_grad(true);    // Not yet implemented
+    //m_dec.has_grad(true);   // Not yet implemented
+    m_radius.has_grad(true);
  
     // Set parameter pointer(s)
     m_pars.push_back(&m_radius);
  
     // Initialise precomputation cache. Note that zero values flag
     // uninitialised as a zero radius is not meaningful
-    m_last_radius = 0.0;
-    m_radius_rad  = 0.0;
-    m_norm        = 0.0;
+    m_last_radius     = 0.0;
+    m_radius_rad      = 0.0;
+    m_norm            = 0.0;
+    m_factor_gradient = 0.0;
  
     // Return
     return;
@@ -554,9 +585,10 @@ void GModelSpatialRadialDisk::copy_members(const GModelSpatialRadialDisk&amp; model)
     m_radius = model.m_radius;
  
     // Copy precomputation cache
-    m_last_radius = model.m_last_radius;
-    m_radius_rad  = model.m_radius_rad;
-    m_norm        = model.m_norm;
+    m_last_radius     = model.m_last_radius;
+    m_radius_rad      = model.m_radius_rad;
+    m_norm            = model.m_norm;
+    m_factor_gradient = model.m_factor_gradient;
  
     // Return
     return;
@@ -576,6 +608,8 @@ void GModelSpatialRadialDisk::free_members(void)
 /***********************************************************************//**
  * @brief Update precomputation cache
  *
+ * @param[in] gradients Request gradient computation?
+ *
  * Computes the normalization
  * \f[
  *    {\tt m\_norm} = \frac{1}{2 \pi (1 - \cos r)}
@@ -586,7 +620,7 @@ void GModelSpatialRadialDisk::free_members(void)
  * approximation you might have expected:
  * \f[{\tt m\_norm} = \frac{1}{\pi r ^ 2}\f]
  ***************************************************************************/
-void GModelSpatialRadialDisk::update() const
+void GModelSpatialRadialDisk::update(const bool& gradients) const
 {
     // Update if radius has changed
     if (m_last_radius != radius()) {
@@ -597,10 +631,26 @@ void GModelSpatialRadialDisk::update() const
         // Compute disk radius in radians
         m_radius_rad = radius() * gammalib::deg2rad;
  
-        // Perform precomputations
-        double denom = gammalib::twopi * (1 - std::cos(m_radius_rad));
+        // Compute sine and cosine of radius, benefitting of the sincos
+        // function
+        double sin_radius;
+        double cos_radius;
+        if (gradients) {
+            sin_radius = std::sin(m_radius_rad);
+            cos_radius = std::cos(m_radius_rad);
+        }
+        else {
+            cos_radius = std::cos(m_radius_rad);
+        }
+
+        // Precompute radial disk function
+        double denom = gammalib::twopi * (1.0 - cos_radius);
         m_norm       = (denom > 0.0) ? 1.0 / denom : 0.0;
  
+        // Precompute factor gradient
+        m_factor_gradient = -gammalib::twopi * m_norm * m_norm * sin_radius *
+                            m_radius.scale() * gammalib::deg2rad;
+
     } // endif: update required
  
     // Return
+#! /usr/bin/env python
+# ==========================================================================
+# Approximate radial disk function
+# --------------------------------------------------------------------------
+#
+# Copyright (C) 2020 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import math
+import gammalib
+import matplotlib.pyplot as plt
+
+
+# ========================== #
+# Compare radial model to MC #
+# ========================== #
+def show(theta0, k=1000.0, rad_max=1.0, nrad=1000):
+    """
+    Show smooth approximation of radial disk function.
+    """
+    # Get dummy energy and time
+    energy = gammalib.GEnergy()
+    time   = gammalib.GTime()
+    
+    # Generate theta vector
+    theta_bin = rad_max/nrad
+    theta     = [(i+0.5) * theta_bin for i in range(nrad)]
+
+    # Generate model vector
+    model = [1.0 / (1.0 + math.exp(-k * (theta0 - t))) for t in theta]
+    disk  = [1.0 if t < theta0 else 0.0 for t in theta]
+
+    # Create figure
+    plt.figure(1)
+    plt.title('Smooth approximation of radial disk function')
+
+    # Plot data
+    plt.plot(theta, model, 'r-')
+    plt.plot(theta, disk, 'g-')
+
+    # Set axes
+    plt.xlabel('Offset (deg)')
+    plt.ylabel('Arbitrary units')
+
+    # Notify
+    print('PLEASE CLOSE WINDOW TO CONTINUE ...')
+
+    # Show plot
+    plt.show()
+
+    # Return
+    return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+    # Show radial model
+    show(0.5)