use simple model derivative computation

knodlseder
1 parent a6683935
Showing 2 changed files with 14 additions and 8 deletions
include/GObservation.hpp
src/obs/GObservation.cpp
@@ -90,14 +90,16 @@ public:
     const GEvents*        events(void) const;
     const std::string&    statistics(void) const { return m_statistics; }
  
+    // Other methods
+    virtual double model_grad(const GModel& model, const GEvent& event, int ipar) const;
+    virtual double npred_grad(const GModel& model, int ipar) const;
+
 protected:
     // Protected methods
     void init_members(void);
     void copy_members(const GObservation& obs);
     void free_members(void);
  
-    // Model methods
-    virtual double model_grad(const GModel& model, const GEvent& event, int ipar) const;
  
     // Model gradient kernel classes
     class model_func : public GFunction {
@@ -119,7 +121,6 @@ protected:
     };
  
     // Npred methods
-    virtual double npred_grad(const GModel& model, int ipar) const;
     virtual double npred_temp(const GModel& model) const;
     virtual double npred_spec(const GModel& model, const GTime& obsTime) const;
  
@@ -425,6 +425,12 @@ void GObservation::free_members(void)
  * @param[in] event Event.
  * @param[in] ipar Parameter index for which gradient should be returned.
  *
+ * This method uses a robust but dumb method to estimate parameter
+ * gradients that have not been provided by the model. We use here a dumb
+ * method as this method is likely used for spatial model parameters, and
+ * the spatial model may eventually be noisy due to numerical integration
+ * limits.
+ *
  * @todo We simply remove any parameter boundaries here for the computation
  *       to avoid any out of boundary errors. We may have models, however,
  *       for which out of bound parameters lead to illegal computations, such
@@ -461,13 +467,12 @@ double GObservation::model_grad(const GModel&amp; model, const GEvent&amp; event,
             // Setup derivative function
             GObservation::model_func function(this, model, event, ipar);
  
-            // Get derivative
+            // Get derivative. We use a fixed step size here that has been
+            // checked on spatial parameters of models
             GDerivative derivative(&function);
             double x  = model[ipar].value();
-            double dx = 0.1;
-            derivative.eps(1.0e-4);
-            //derivative.silent(true);
-            grad = derivative.value(x, dx);
+            double dx = 0.05;
+            grad = derivative.difference(x, dx);
  
             // Restore current model parameter
             (*ptr)[ipar] = current;