diff --git a/doc/dev/inst/cta/gammalib_cta.tex b/doc/dev/inst/cta/gammalib_cta.tex index 3f80440..745a3ec 100644 --- a/doc/dev/inst/cta/gammalib_cta.tex +++ b/doc/dev/inst/cta/gammalib_cta.tex @@ -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} diff --git a/include/GModelSpatialRadialDisk.hpp b/include/GModelSpatialRadialDisk.hpp index 6aac673..01a851d 100644 --- a/include/GModelSpatialRadialDisk.hpp +++ b/include/GModelSpatialRadialDisk.hpp @@ -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 }; diff --git a/src/model/GModelSpatialRadialDisk.cpp b/src/model/GModelSpatialRadialDisk.cpp index 304a769..2cd0ce4 100644 --- a/src/model/GModelSpatialRadialDisk.cpp +++ b/src/model/GModelSpatialRadialDisk.cpp @@ -51,6 +51,7 @@ const GModelSpatialRegistry g_radial_disk_legacy_registry(&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& 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& 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 diff --git a/test/dev/test_disk_approximation.py b/test/dev/test_disk_approximation.py new file mode 100755 index 0000000..493cd4e --- /dev/null +++ b/test/dev/test_disk_approximation.py @@ -0,0 +1,74 @@ +#! /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 . +# +# ========================================================================== +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)