Merge branch 'devel' of https://cta-gitlab.irap.omp.eu/gammalib/gammalib into devel

ldivenere
2 parents a54c4b50 bd372e66
Showing 10 changed files with 354 additions and 43 deletions
AUTHORS
ChangeLog
NEWS
doc/source/admin/release_history/1.6.rst
include/GWcs.hpp
inst/cta/src/GCTAOnOffObservation.cpp
pyext/GPythonOptimizerFunction.i
pyext/gammalib/opt.i
src/sky/GSkyDir.cpp
src/sky/GWcs.cpp
 Rolf Buehler (rolf.buehler@desy.de)
-François Brun (brun@cenbg.in2p3.fr)
+François Brun (francois.brun@cea.fr)
 Joshua Cardenzana (jcardenzana@irap.omp.eu)
 Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr)
 Christoph Deil (Christoph.Deil@mpi-hd.mpg.de)
-2018-11-27
+2019-01-09
  
         * Version 1.6.0 released
         ========================
  
+        Improve OpenMP handling in GWcs::dir2pix (#2770)
         Add support for Fermi/LAT P8V3 response functions
         CTA background rates are now per ontime (#2758)
+        Add GPythonOptimizerFunction class (#2733)
         Remove diffuse response cube caching (#2722)
         Revise GCTAEdispCube::operator() interface and disable stacked response
           cutting in true energy (#2721)
 New Features and Important Changes in GammaLib 1.6.0
  
-23 November 2018
+9 January 2019
  
  
 1. Introduction
@@ -14,6 +14,7 @@ made since the last release of GammaLib.
 The following classes have been added:
 - GModelSpectralExponential
 - GWcsARC
+- GPythonOptimizerFunction
 - GCTAResponseCache
 - GCTAModelBackground
 - GCTAModelSpatial
@@ -247,6 +248,9 @@ None
  
 12. Sky module
 --------------
+An OpenMP locking mechanism was added to GWcs::dir2pix that improves the speed of
+the code parallelization by avoiding mutual lockups of the threads (#2770).
+
 The class GWcsARC was added to handle ARC WCS projections (#2541). In addition
 the private GWcs::prj_bchk() method was added to implement boundary checking
 similar to the checking done in WCS.
@@ -49,7 +49,11 @@ Bug fixes
 Improvements
 ------------
  
+* [`2770 <https://cta-redmine.irap.omp.eu/issues/2770>`_] -
+  Improve OpenMP handling in ``GWcs::dir2pix``
 * Add support for Fermi/LAT P8V3 response functions
+* [`2733 <https://cta-redmine.irap.omp.eu/issues/2733>`_] -
+  Add ``GPythonOptimizerFunction`` class
 * [`2721 <https://cta-redmine.irap.omp.eu/issues/2721>`_] -
   Revise ``GCTAEdispCube::operator()`` interface and disable stacked response
   cutting in true energy
 /***************************************************************************
  *          GWcs.hpp - Abstract world coordinate system base class         *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2018 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2019 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -30,6 +30,7 @@
 /* __ Includes ___________________________________________________________ */
 #include <vector>
 #include <string>
+#include <map>
 #include <iostream>
 #include "GSkyProjection.hpp"
 #include "GLog.hpp"
@@ -37,6 +38,10 @@
 #include "GSkyDir.hpp"
 #include "GSkyPixel.hpp"
  
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
  
 /***********************************************************************//**
  * @class GWcs
@@ -220,6 +225,15 @@ protected:
     mutable GSkyDir   m_last_dir2pix_dir;  //!< Last sky direction for dir2pix
     mutable GSkyPixel m_last_pix2dir_pix;  //!< Last pixel for pix2dir
     mutable GSkyPixel m_last_dir2pix_pix;  //!< Last pixel for dir2pix
+    
+    // Thread locking
+    void init_lock(const int& lock_id=0) const;
+    void set_lock(const int& lock_id=0) const;
+    void unset_lock(const int& lock_id=0) const;
+
+    #ifdef _OPENMP
+    mutable std::map<int, omp_lock_t> m_locks;
+    #endif
 };
  
  
@@ -236,4 +250,49 @@ int GWcs::size(void) const
     return 2;
 }
  
+
+/***********************************************************************//**
+ * @brief Initializes an OpenMP lock with a specific name
+ *
+ * @param[in] lock_id Identifier of the lock that should be initialized
+ ***************************************************************************/
+inline
+void GWcs::init_lock(const int& lock_id) const
+{
+    #ifdef _OPENMP
+    if (m_locks.count(lock_id) == 0) {
+        m_locks[lock_id] = omp_lock_t();
+        omp_init_lock(&m_locks[lock_id]);
+    }
+    #endif
+}
+
+
+/***********************************************************************//**
+ * @brief Sets an OpenMP lock with a specific name
+ *
+ * @param[in] lock_id Identifier of the lock that should be set
+ ***************************************************************************/
+inline
+void GWcs::set_lock(const int& lock_id) const
+{
+    #ifdef _OPENMP
+    omp_set_lock(&m_locks[lock_id]);
+    #endif
+}
+
+
+/***********************************************************************//**
+ * @brief Releases a previously set OpenMP lock 
+ * 
+ * @param[in] lock_id Identifier of the lock to be released
+ ***************************************************************************/
+inline
+void GWcs::unset_lock(const int& lock_id) const
+{
+    #ifdef _OPENMP
+    omp_unset_lock(&m_locks[lock_id]);
+    #endif
+}
+
 #endif /* GWCS_HPP */
@@ -2675,8 +2675,11 @@ double GCTAOnOffObservation::likelihood_cstat(const GModels&amp; models,
         double ngam = N_gamma(models, i, &sky_grad);
         double nbgd = N_bgd(models, i, &bgd_grad);
  
-        // Skip bin if model is too small (avoids -Inf or NaN gradients)
+        // Compute and updated predicted number of events
         double nonpred = ngam + alpha * nbgd;
+        *npred        += nonpred;
+
+        // Skip bin if model is too small (avoids -Inf or NaN gradients)
         if ((nbgd <= minmod) || (nonpred <= minmod)) {
             #if defined(G_LIKELIHOOD_DEBUG)
             n_small_model++;
@@ -2685,10 +2688,8 @@ double GCTAOnOffObservation::likelihood_cstat(const GModels&amp; models,
         }
  
         // Now we have all predicted gamma and background events for
-        // current energy bin. Update the log(likelihood) and predicted
-        // number of events
-        value  += -non * log(nonpred) + nonpred - noff * log(nbgd) + nbgd;
-        *npred += nonpred;
+        // current energy bin. Update the log(likelihood)
+        value += -non * std::log(nonpred) + nonpred - noff * std::log(nbgd) + nbgd;
  
         // Update statistics
         #if defined(G_LIKELIHOOD_DEBUG)
+/***************************************************************************
+ *       GPythonOptimizerFunction.i - Python optimizer function class      *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2018 by 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/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file GPythonOptimizerFunction.i
+ * @brief Python optimizer function class
+ * @author Juergen Knoedlseder
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "GOptimizerFunction.hpp"
+#include "GTools.hpp"
+#include "GVector.hpp"
+#include "GMatrixSparse.hpp"
+#include "GException.hpp"
+
+
+/***********************************************************************//**
+ * @class GPythonException
+ *
+ * @brief Python exception
+ *
+ * This class implements a Python exception. It allows to catch any error
+ * message that occurs during function evaluation.
+ ***************************************************************************/
+class GPythonException : public GExceptionHandler {
+public:
+    class eval_error : public GExceptionHandler {
+    public:
+        eval_error(const std::string& origin) {
+            
+            // Set origin and default message
+            m_origin  = origin;
+            m_message = "";
+
+            // Set message
+            if (PyErr_Occurred()) {
+
+                // Fetch error type, value and traceback
+                PyObject *type      = 0;
+                PyObject *value     = 0;
+                PyObject *traceback = 0;
+                PyErr_Fetch(&type, &value, &traceback);
+
+                // Clear error message
+                PyErr_Clear();
+
+                // Extract message strings
+                PyObject *py_type  = PyObject_Str(type);
+                PyObject *py_value = PyObject_Str(value);
+                char     *c_type   = NULL;
+                char     *c_value  = NULL;
+                if (PyUnicode_Check(py_type)) {
+                    PyObject* temp_bytes = PyUnicode_AsEncodedString(py_type,
+                                                                     "utf-8",
+                                                                     "Error ~");
+                    if (temp_bytes != NULL) {
+                        c_type = PyString_AsString(temp_bytes);
+                    }
+                }
+                else if (py_type != NULL) {
+                    c_type = PyString_AsString(py_type);
+                }
+                if (c_type != NULL) {
+                    m_message += std::string(c_type);
+                    m_message += "\n";
+                }
+                if (PyUnicode_Check(py_value)) {
+                    PyObject* temp_bytes = PyUnicode_AsEncodedString(py_value,
+                                                                     "utf-8",
+                                                                     "Error ~");
+                    if (temp_bytes != NULL) {
+                        c_value = PyString_AsString(temp_bytes);
+                    }
+                }
+                else if (py_value != NULL) {
+                    c_value = PyString_AsString(py_value);
+                }
+                if (c_value != NULL) {
+                    m_message += std::string(c_value);
+                }
+                Py_DECREF(py_type);
+                Py_DECREF(py_value);
+            }
+        }
+    };
+};
+
+
+/***********************************************************************//**
+ * @class GPythonOptimizerFunction
+ *
+ * @brief Python interface for optimizer function
+ *
+ * This class implements the call back for a Python optimizer function,
+ * allowing to use the GammaLib optimizer classes from within Python.
+ ***************************************************************************/
+class GPythonOptimizerFunction : public GOptimizerFunction {
+public:
+    // Constructor
+    GPythonOptimizerFunction(void) {
+        m_py_object = NULL;
+        m_value     = 0.0;
+        m_gradient.clear();
+        m_curvature.clear();
+        m_pars.clear();
+    }
+
+    // Destructor
+    virtual ~GPythonOptimizerFunction(void) {
+        if (m_py_object != NULL) {
+            Py_DECREF(m_py_object);
+        }
+    }
+
+    // Evaluate Python function
+    void eval(const GOptimizerPars& pars) {
+
+        // Store parameters
+        m_pars = pars;
+
+        // Initialise value, gradient and curvature
+        int npars = pars.size();
+        if (npars < 1) {
+            npars = 1;
+        }
+
+        // Initialise value, gradient vector and curvature matrix
+        m_value     = 0.0;
+        m_gradient  = GVector(npars);
+        m_curvature = GMatrixSparse(npars,npars);
+
+        // Execute evaluation
+        PyObject* args = Py_BuildValue("()");
+        PyObject* func = static_cast<PyObject*>(m_py_object);
+        PyObject* res  = PyEval_CallObject(func, args);
+        Py_DECREF(args);
+
+        // Release result
+        if (res != NULL) {
+            Py_DECREF(res);
+        }
+        else {
+            throw GPythonException::eval_error("GPythonOptimizerFunction::eval");
+        }
+    }
+
+    // Return value
+    double value(void) const {
+        return m_value;
+    }
+
+    // Return gradient
+    GVector* gradient(void) {
+        return &m_gradient;
+    }
+
+    // Return curvature
+    GMatrixSparse* curvature(void) {
+        return &m_curvature;
+    }
+
+    // Return parameters
+    const GOptimizerPars& _pars(void) const {
+        return m_pars;
+    }
+
+    // Set Python callback function
+    void _set_eval(PyObject* function) {
+        m_py_object = function;
+        Py_INCREF(function);
+    }
+
+    // Set value
+    void _set_value(const double& value) {
+        m_value = value;
+    }
+
+    // Callback information
+    PyObject*      m_py_object;  //!< Python callback function
+    double         m_value;      //!< Function value
+    GVector        m_gradient;   //!< Function gradient
+    GMatrixSparse  m_curvature;  //!< Function curvature
+    GOptimizerPars m_pars;       //!< Current parameters
+};
+%}
+
+
+/***********************************************************************//**
+ * @class GPythonOptimizerFunction
+ *
+ * @brief Optimizer class Python interface definition
+ ***************************************************************************/
+class GPythonOptimizerFunction : public GOptimizerFunction {
+public:        
+    // Constructors and destructors
+    GPythonOptimizerFunction(void);
+    virtual ~GPythonOptimizerFunction(void);
+
+    // Methods
+    void           eval(const GOptimizerPars& pars);
+    double         value(void) const;
+    GVector*       gradient(void);
+    GMatrixSparse* curvature(void);
+
+    // Protected methods
+    const GOptimizerPars& _pars(void) const;
+    void                  _set_eval(PyObject* function);
+    void                  _set_value(const double& value);
+};
 /***************************************************************************
  *                       opt.i - Optimization module                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -90,3 +90,4 @@
 %include "GOptimizerPar.i"
 %include "GOptimizerPars.i"
 %include "GOptimizerFunction.i"
+%include "GPythonOptimizerFunction.i"
@@ -800,51 +800,52 @@ bool operator==(const GSkyDir &amp;a, const GSkyDir &amp;b)
  
     // Compute dependent on coordinate system availability. This speeds
     // up things by avoiding unnecessary coordinate transformations.
-    if (a.m_has_lb) {
-        if (b.m_has_lb) {
-            if (std::abs(a.m_b) == 90.0) {
-                equal = (a.m_b == b.m_b);
-            }
-            else {
-                equal = (a.m_b == b.m_b && a.m_l == b.m_l);
-            }
+
+    // Check if both have equatorial coordinates
+    if (a.m_has_radec && b.m_has_radec) {
+        if (std::abs(a.m_dec) == 90.0) {
+            equal = (a.m_dec == b.m_dec);
         }
         else {
+            equal = (a.m_dec == b.m_dec && a.m_ra == b.m_ra);
+        }
+    }
+    // ... check if both have Galactic coordinates
+    else if (a.m_has_lb && b.m_has_lb) {
+        if (std::abs(a.m_b) == 90.0) {
+            equal = (a.m_b == b.m_b);
+        }
+        else {
+            equal = (a.m_b == b.m_b && a.m_l == b.m_l);
+        }
+    }
+    // ... otherwise the coordinate systems are different
+    else {
+        if (a.m_has_lb) {
             if (std::abs(a.m_b) == 90.0) {
                 equal = (a.m_b == b.b());
             }
             else {
                 equal = (a.m_b == b.b() && a.m_l == b.l());
             }
-        }
-    }
-    else if (a.m_has_radec) {
-        if (b.m_has_radec) {
+        } 
+        else if (a.m_has_radec) {
             if (std::abs(a.m_dec) == 90.0) {
-                equal = (a.m_dec == b.m_dec);
+                equal = (a.m_dec == b.dec());
             }
             else {
-            equal = (a.m_dec == b.m_dec && a.m_ra == b.m_ra);
+                equal = (a.m_dec == b.dec() && a.m_ra == b.ra());
             }
         }
         else {
-            if (std::abs(a.m_dec) == 90.0) {
-            equal = (a.m_dec == b.dec());
+            if (std::abs(b.dec()) == 90.0) {
+                equal = (b.dec() == a.dec());
             }
             else {
-                equal = (a.m_dec == b.dec() && a.m_ra == b.ra());
+                equal = (b.dec() == a.dec() && b.ra() == a.ra());
             }
         }
     }
-    else {
-        if (std::abs(b.dec()) == 90.0) {
-            equal = (b.dec() == a.dec());
-        }
-        else {
-            equal = (b.dec() == a.dec() && b.ra() == a.ra());
-        }
-    }
-
     // Return equality
     return equal;
 }
 /***************************************************************************
  *          GWcs.cpp - Abstract world coordinate system base class         *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2018 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2019 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -678,13 +678,16 @@ GSkyDir GWcs::pix2dir(const GSkyPixel&amp; pixel) const
     // the one used for the cache, then get the cached sky direction. We need
     // to put this in a OMP critical zone so that no other thread will change
     // the cache in the meantime
-    #pragma omp critical(GWcs_pix2dir_1)
+    bool update_cache = true;
+    set_lock(2);
     if (m_has_pix2dir_cache && (m_last_pix2dir_pix == pixel)) {
         dir = m_last_pix2dir_dir;
+        update_cache = false;
     }
+    unset_lock(2);
  
     // ... otherwise the computation of the sky direction is needed
-    else {
+    if (update_cache) {
  
         // Allocate memory for transformation
         double pixcrd[2];
@@ -712,12 +715,13 @@ GSkyDir GWcs::pix2dir(const GSkyPixel&amp; pixel) const
  
         // Store result in cache. We need to put this in a OMP critical zone
         // to make sure that no other thread is fiddeling with the cache
-        #pragma omp critical(GWcs_pix2dir_2)
+        set_lock(2);
         {
             m_has_pix2dir_cache = true;
             m_last_pix2dir_pix  = pixel;
             m_last_pix2dir_dir  = dir;
         }
+        unset_lock(2);
  
         // Debug: Dump transformation steps
         #if defined(G_XY2DIR_DEBUG)
@@ -758,13 +762,16 @@ GSkyPixel GWcs::dir2pix(const GSkyDir&amp; dir) const
     // to the one used for the cache, then get the cached sky pixel. We need
     // to put this in a OMP critical zone so that no other thread will change
     // the cache in the meantime
-    #pragma omp critical(GWcs_dir2pix_1)
+    bool update_cache = true;
+    set_lock(1);
     if (m_has_dir2pix_cache && (m_last_dir2pix_dir == dir)) {
         pixel = m_last_dir2pix_pix;
+        update_cache = false;
     }
+    unset_lock(1);
  
     // ... otherwise the computation of the sky pixel is needed
-    else {
+    if (update_cache) {
  
         // Allocate memory for transformation
         double pixcrd[2];
@@ -793,12 +800,13 @@ GSkyPixel GWcs::dir2pix(const GSkyDir&amp; dir) const
  
         // Store result in cache. We need to put this in a OMP critical zone
         // to make sure that no other thread is fiddeling with the cache
-        #pragma omp critical(GWcs_dir2pix_2)
+        set_lock(1);
         {
             m_has_dir2pix_cache = true;
             m_last_dir2pix_dir  = dir;
             m_last_dir2pix_pix  = pixel;
         }
+        unset_lock(1);
  
         // Debug: Dump transformation steps
         #if defined(G_DIR2XY_DEBUG)
@@ -944,6 +952,10 @@ void GWcs::init_members(void)
     m_last_pix2dir_pix.clear();
     m_last_dir2pix_pix.clear();
  
+    // Initialize the locks
+    init_lock(1);
+    init_lock(2);
+
     // Return
     return;
 }