From c135e7ceca61f203c46c583ce1d1bc3e31ad7ffa Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 26 Sep 2016 14:28:31 +0200 Subject: [PATCH] Remove eval_gradients() methods and remplace by eval() method with optional gradients argument (#1847) --- ChangeLog | 1 + NEWS | 10 ++++++++-- doc/source/users/user_manual/modules/model.rst | 4 ++-- doc/source/users/user_manual/modules/obs.rst | 6 +++--- include/GModel.hpp | 16 +++++++--------- include/GModelData.hpp | 7 +++---- include/GModelSky.hpp | 11 ++++------- include/GModelSpatial.hpp | 4 ++-- include/GModelSpatialDiffuse.hpp | 4 ++-- include/GModelSpatialDiffuseConst.hpp | 4 ++-- include/GModelSpatialDiffuseCube.hpp | 4 ++-- include/GModelSpatialDiffuseMap.hpp | 4 ++-- include/GModelSpatialElliptical.hpp | 15 ++++++--------- include/GModelSpatialEllipticalDisk.hpp | 7 ++----- include/GModelSpatialEllipticalGauss.hpp | 7 ++----- include/GModelSpatialPointSource.hpp | 4 ++-- include/GModelSpatialRadial.hpp | 14 ++++++-------- include/GModelSpatialRadialDisk.hpp | 6 ++---- include/GModelSpatialRadialGauss.hpp | 6 ++---- include/GModelSpatialRadialShell.hpp | 6 ++---- include/GModelSpectral.hpp | 5 ++--- include/GModelSpectralBrokenPlaw.hpp | 5 ++--- include/GModelSpectralConst.hpp | 5 ++--- include/GModelSpectralExpInvPlaw.hpp | 5 ++--- include/GModelSpectralExpPlaw.hpp | 5 ++--- include/GModelSpectralFunc.hpp | 5 ++--- include/GModelSpectralGauss.hpp | 5 ++--- include/GModelSpectralLogParabola.hpp | 5 ++--- include/GModelSpectralNodes.hpp | 5 ++--- include/GModelSpectralPlaw.hpp | 5 ++--- include/GModelSpectralPlawEnergyFlux.hpp | 5 ++--- include/GModelSpectralPlawPhotonFlux.hpp | 5 ++--- include/GModelSpectralSuperExpPlaw.hpp | 5 ++--- include/GModelTemporal.hpp | 6 +++--- include/GModelTemporalConst.hpp | 6 +++--- include/GModels.hpp | 17 +++++++++++------ include/GOptimizerPar.hpp | 30 +++++++++++++++--------------- inst/com/include/GCOMModelDRBFitting.hpp | 7 +++---- inst/com/pyext/GCOMModelDRBFitting.i | 7 +++---- inst/com/src/GCOMModelDRBFitting.cpp | 161 ++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------- inst/cta/include/GCTAModelAeffBackground.hpp | 7 +++---- inst/cta/include/GCTAModelCubeBackground.hpp | 7 +++---- inst/cta/include/GCTAModelIrfBackground.hpp | 7 +++---- inst/cta/include/GCTAModelRadial.hpp | 6 +++--- inst/cta/include/GCTAModelRadialAcceptance.hpp | 7 +++---- inst/cta/include/GCTAModelRadialGauss.hpp | 35 ++++++++++++++++++++++++++++++----- inst/cta/include/GCTAModelRadialPolynom.hpp | 6 +++--- inst/cta/include/GCTAModelRadialProfile.hpp | 6 +++--- inst/cta/pyext/GCTAModelAeffBackground.i | 31 +++++++++++++++---------------- inst/cta/pyext/GCTAModelCubeBackground.i | 7 +++---- inst/cta/pyext/GCTAModelIrfBackground.i | 7 +++---- inst/cta/pyext/GCTAModelRadial.i | 6 +++--- inst/cta/pyext/GCTAModelRadialAcceptance.i | 7 +++---- inst/cta/pyext/GCTAModelRadialGauss.i | 6 +++--- inst/cta/pyext/GCTAModelRadialPolynom.i | 6 +++--- inst/cta/pyext/GCTAModelRadialProfile.i | 6 +++--- inst/cta/src/GCTAModelAeffBackground.cpp | 125 +++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------ inst/cta/src/GCTAModelCubeBackground.cpp | 106 +++++++++++++++++++++++++++++----------------------------------------------------------------------------- inst/cta/src/GCTAModelIrfBackground.cpp | 120 ++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------ inst/cta/src/GCTAModelRadialAcceptance.cpp | 138 ++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------ inst/cta/src/GCTAModelRadialGauss.cpp | 66 +++++++++++++++++++----------------------------------------------- inst/cta/src/GCTAModelRadialPolynom.cpp | 96 ++++++++++++++++++++++++++---------------------------------------------------------------------- inst/cta/src/GCTAModelRadialProfile.cpp | 27 +++++++-------------------- inst/cta/src/GCTAOnOffObservation.cpp | 16 ++++++++++------ inst/cta/test/dev/test_gauss.py | 4 ++-- inst/cta/test/dev/test_radial_acceptance.py | 2 +- pyext/GModel.i | 7 +++---- pyext/GModelData.i | 7 +++---- pyext/GModelSky.i | 5 ++--- pyext/GModelSpatial.i | 6 +++--- pyext/GModelSpatialDiffuse.i | 6 +++--- pyext/GModelSpatialDiffuseConst.i | 6 +++--- pyext/GModelSpatialDiffuseCube.i | 8 ++++---- pyext/GModelSpatialDiffuseMap.i | 4 ++-- pyext/GModelSpatialElliptical.i | 17 +++++++---------- pyext/GModelSpatialEllipticalDisk.i | 16 +++++----------- pyext/GModelSpatialEllipticalGauss.i | 16 +++++----------- pyext/GModelSpatialPointSource.i | 6 +++--- pyext/GModelSpatialRadial.i | 12 +++++------- pyext/GModelSpatialRadialDisk.i | 21 ++++++++------------- pyext/GModelSpatialRadialGauss.i | 21 ++++++++------------- pyext/GModelSpatialRadialShell.i | 21 ++++++++------------- pyext/GModelSpectral.i | 5 ++--- pyext/GModelSpectralBrokenPlaw.i | 5 ++--- pyext/GModelSpectralConst.i | 5 ++--- pyext/GModelSpectralExpInvPlaw.i | 8 ++------ pyext/GModelSpectralExpPlaw.i | 5 ++--- pyext/GModelSpectralFunc.i | 5 ++--- pyext/GModelSpectralGauss.i | 7 +++---- pyext/GModelSpectralLogParabola.i | 5 ++--- pyext/GModelSpectralNodes.i | 5 ++--- pyext/GModelSpectralPlaw.i | 5 ++--- pyext/GModelSpectralPlawEnergyFlux.i | 5 ++--- pyext/GModelSpectralPlawPhotonFlux.i | 5 ++--- pyext/GModelSpectralSuperExpPlaw.i | 5 ++--- pyext/GModelTemporal.i | 12 +++++++----- pyext/GModelTemporalConst.i | 12 +++++++----- pyext/GModels.i | 5 +++-- src/model/GModelSky.cpp | 43 ++++++++++++------------------------------- src/model/GModelSpatialDiffuseConst.cpp | 55 ++++++++++++++++--------------------------------------- src/model/GModelSpatialDiffuseCube.cpp | 52 +++++++++++++++++++--------------------------------- src/model/GModelSpatialDiffuseMap.cpp | 37 +++++++++++++------------------------ src/model/GModelSpatialElliptical.cpp | 40 +++++++++++----------------------------- src/model/GModelSpatialEllipticalDisk.cpp | 33 +++++++-------------------------- src/model/GModelSpatialEllipticalGauss.cpp | 33 +++++++-------------------------- src/model/GModelSpatialPointSource.cpp | 29 +++++++---------------------- src/model/GModelSpatialRadial.cpp | 33 +++++++++------------------------ src/model/GModelSpatialRadialDisk.cpp | 31 +++++++------------------------ src/model/GModelSpatialRadialGauss.cpp | 28 +++++++--------------------- src/model/GModelSpatialRadialShell.cpp | 30 ++++++++---------------------- src/model/GModelSpectralBrokenPlaw.cpp | 134 +++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------- src/model/GModelSpectralConst.cpp | 52 +++++++++++++++------------------------------------- src/model/GModelSpectralExpInvPlaw.cpp | 113 +++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------ src/model/GModelSpectralExpPlaw.cpp | 113 +++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------ src/model/GModelSpectralFunc.cpp | 72 ++++++++++++++++-------------------------------------------------------- src/model/GModelSpectralGauss.cpp | 97 ++++++++++++++++++++++++++----------------------------------------------------------------------- src/model/GModelSpectralLogParabola.cpp | 108 ++++++++++++++++++++++++++++++------------------------------------------------------------------------------ src/model/GModelSpectralNodes.cpp | 103 ++++++++++++++++++++++++++----------------------------------------------------------------------------- src/model/GModelSpectralPlaw.cpp | 101 ++++++++++++++++++++++++++++------------------------------------------------------------------------- src/model/GModelSpectralPlawEnergyFlux.cpp | 107 ++++++++++++++++++++++------------------------------------------------------------------------------------- src/model/GModelSpectralPlawPhotonFlux.cpp | 108 +++++++++++++++++++++--------------------------------------------------------------------------------------- src/model/GModelSpectralSuperExpPlaw.cpp | 129 +++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------------------- src/model/GModelTemporalConst.cpp | 49 ++++++++++++++++--------------------------------- src/model/GModels.cpp | 36 +++++++++--------------------------- src/obs/GObservation.cpp | 7 +------ src/obs/GResponse.cpp | 69 +++++++++++++++++++++++---------------------------------------------- test/testinst/GTestModelData.hpp | 14 +++++--------- 127 files changed, 1079 insertions(+), 2351 deletions(-) diff --git a/ChangeLog b/ChangeLog index 17b525c..efbda85 100644 --- a/ChangeLog +++ b/ChangeLog @@ -4,6 +4,7 @@ ========================= Check CTA response table dimension for consistency (#1848) + Remove eval_gradients() methods and remplace by eval() method with optional gradients argument (#1847) Add GTestSuite::test_value() method with automatic precision computation Add GSkyMap::is_empty() method (#1828) Add GApplication access by index operator (#1838) diff --git a/NEWS b/NEWS index 3450da8..f98c080 100644 --- a/NEWS +++ b/NEWS @@ -36,14 +36,20 @@ The following methods have been added: - GSkyMap::is_empty() The following methods have been removed: -- ... +- ::eval_gradients() + The eval_gradients() methods duplicate much of the code of the eval() + methods. To increase the maintainability of the code, the eval_gradients() + methods have been removed, and the computation of partial derivatives is + now done by the eval() methods if the optional gradients argument is set + to true. The following methods have been renamed: - GModelSpectralPlawPhotonFlux::flux(), was GModelSpectralPlaw2::integral() - GLog::buffer_size(), was GLog::max_size() -The parameters for the following methods have been changed: +The arguments for the following methods have been changed: - GApplication::log_parameters() now has a mandatory GChatter argument +- ::eval() has now an optional gradients parameter (see above) The return value of the following methods has been changed: - GSkyMap::write() now returns a pointer to GFitsHDU instead of void diff --git a/doc/source/users/user_manual/modules/model.rst b/doc/source/users/user_manual/modules/model.rst index c28db07..0e101c7 100644 --- a/doc/source/users/user_manual/modules/model.rst +++ b/doc/source/users/user_manual/modules/model.rst @@ -35,7 +35,7 @@ methods. For model evaluation, the :doxy:``GObservation::model`` method loops over all models in the :doxy:`GModels` container and calls the -:doxy:`GModel::eval_gradients` method of each model. +:doxy:`GModel::eval` method of each model. For :doxy:`GModelSky` models this method invokes the :doxy:`GResponse::convolve` method that calls among others the :doxy:`GResponse::irf` method. That method needs to be implement for each instruments. @@ -46,7 +46,7 @@ The integrate over observed time is for the moment a simple multiplication with ontime. The integral over observed energy is performed by the ``GObservation::npred_spec`` method that calls the -:doxy:`GModel::eval_gradients` method. +:doxy:`GModel::eval` method. For :doxy:`GModelSky` models this method calls the :doxy:`GResponse::nroi` method that needs to be implement for each instruments. diff --git a/doc/source/users/user_manual/modules/obs.rst b/doc/source/users/user_manual/modules/obs.rst index d7ebe74..9e8b9b2 100644 --- a/doc/source/users/user_manual/modules/obs.rst +++ b/doc/source/users/user_manual/modules/obs.rst @@ -233,9 +233,9 @@ takes the :doxy:`GSource` class instead of the :doxy:`GPhoton` class as argument :doxy:`GSource` differs from :doxy:`GPhoton` in that the photon arrival direction :math:`p` is replaced by the spatial component :doxy:`GModelSpatial` of a source model. -Equation :eq:`pirf` is used by the :doxy:`GModelSky::eval` and -:doxy:`GModelSky::eval_gradients` methods for computation of the instrument -response to a source model (see :ref:`fig_calltree_model`). +Equation :eq:`pirf` is used by the :doxy:`GModelSky::eval` method for +computation of the instrument response to a source model +(see :ref:`fig_calltree_model`). A maximum likelihood analysis of the data generally needs the computation of the predicted number of events within the selection region for each source model. diff --git a/include/GModel.hpp b/include/GModel.hpp index 2cfeb46..ee64ca3 100644 --- a/include/GModel.hpp +++ b/include/GModel.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModel.hpp - Abstract virtual model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2015 by Juergen Knoedlseder * + * copyright (C) 2009-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -46,11 +46,10 @@ class GObservation; * * @brief Abstract model class * - * This class implements a parametric model. The class has two methods for - * model evaluation: eval() and eval_gradients(). - * The eval() method evaluates the model for a given event and observation. - * In addition, eval_gradients() also sets the parameter gradients of the - * model. + * This class implements a parametric model. The eval() method evaluates the + * model for a given event and observation. If the gradients parameter is set + * to true, the eval() method also computes analytical parameter gradients if + * they exist. * * A model has the following attributes: * - @p name @@ -118,9 +117,8 @@ public: virtual std::string type(void) const = 0; virtual bool is_constant(void) const = 0; virtual double eval(const GEvent& event, - const GObservation& obs) const = 0; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const = 0; + const GObservation& obs, + const bool& gradients = false) const = 0; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const = 0; virtual void read(const GXmlElement& xml) = 0; diff --git a/include/GModelData.hpp b/include/GModelData.hpp index 951803e..8344f14 100644 --- a/include/GModelData.hpp +++ b/include/GModelData.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelData.hpp - Abstract virtual data model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2013 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -69,9 +69,8 @@ public: virtual std::string type(void) const = 0; virtual bool is_constant(void) const = 0; virtual double eval(const GEvent& event, - const GObservation& obs) const = 0; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const = 0; + const GObservation& obs, + const bool& gradients = false) const = 0; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const = 0; virtual GEvents* mc(const GObservation& obs, GRan& ran) const = 0; diff --git a/include/GModelSky.hpp b/include/GModelSky.hpp index 6f64186..4c35f96 100644 --- a/include/GModelSky.hpp +++ b/include/GModelSky.hpp @@ -92,10 +92,8 @@ class GObservation; * * The class has two methods for model evaluation that evaluate the model * for a specific event, given an observation. The eval() method returns - * the model value, the eval_gradients() returns the model value and sets - * the analytical gradients for all model parameters. - * Note that the eval() and eval_gradients() methods call protected - * methods that handle time dispersion, energy dispersion and the point + * the model value. Note that the eval() method calls protected methods that + * handle time dispersion, energy dispersion and the point * spread function (spatial dispersion). Dispersion is handled by * integrating over the relevant intervals of the source properties. * Integration of the model is done by three nested integrals. @@ -149,9 +147,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; diff --git a/include/GModelSpatial.hpp b/include/GModelSpatial.hpp index 6b341e9..a035d72 100644 --- a/include/GModelSpatial.hpp +++ b/include/GModelSpatial.hpp @@ -87,8 +87,8 @@ public: virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; virtual GClassCode code(void) const = 0; - virtual double eval(const GPhoton& photon) const = 0; - virtual double eval_gradients(const GPhoton& photon) const = 0; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const = 0; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const = 0; virtual double mc_norm(const GSkyDir& dir, diff --git a/include/GModelSpatialDiffuse.hpp b/include/GModelSpatialDiffuse.hpp index 7f68bd3..3147000 100644 --- a/include/GModelSpatialDiffuse.hpp +++ b/include/GModelSpatialDiffuse.hpp @@ -62,8 +62,8 @@ public: virtual GModelSpatialDiffuse* clone(void) const = 0; virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; - virtual double eval(const GPhoton& photon) const = 0; - virtual double eval_gradients(const GPhoton& photon) const = 0; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const = 0; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const = 0; virtual double mc_norm(const GSkyDir& dir, diff --git a/include/GModelSpatialDiffuseConst.hpp b/include/GModelSpatialDiffuseConst.hpp index c5cd85a..d850bb9 100644 --- a/include/GModelSpatialDiffuseConst.hpp +++ b/include/GModelSpatialDiffuseConst.hpp @@ -63,8 +63,8 @@ public: virtual GModelSpatialDiffuseConst* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/include/GModelSpatialDiffuseCube.hpp b/include/GModelSpatialDiffuseCube.hpp index d8585c1..dac9c0a 100644 --- a/include/GModelSpatialDiffuseCube.hpp +++ b/include/GModelSpatialDiffuseCube.hpp @@ -78,8 +78,8 @@ public: virtual GModelSpatialDiffuseCube* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/include/GModelSpatialDiffuseMap.hpp b/include/GModelSpatialDiffuseMap.hpp index 0f3649a..db8795b 100644 --- a/include/GModelSpatialDiffuseMap.hpp +++ b/include/GModelSpatialDiffuseMap.hpp @@ -70,8 +70,8 @@ public: virtual GModelSpatialDiffuseMap* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/include/GModelSpatialElliptical.hpp b/include/GModelSpatialElliptical.hpp index 8b7339a..46868e5 100644 --- a/include/GModelSpatialElliptical.hpp +++ b/include/GModelSpatialElliptical.hpp @@ -69,14 +69,11 @@ public: virtual double eval(const double& theta, const double& posangle, const GEnergy& energy, - const GTime& time) const = 0; - virtual double eval_gradients(const double& theta, - const double& posangle, - const GEnergy& energy, - const GTime& time) const = 0; + const GTime& time, + const bool& gradients = false) const = 0; virtual GSkyDir mc(const GEnergy& energy, - const GTime& time, - GRan& ran) const = 0; + const GTime& time, + GRan& ran) const = 0; virtual bool contains(const GSkyDir& dir, const double& margin = 0.0) const = 0; virtual double theta_max(void) const = 0; @@ -84,8 +81,8 @@ public: // Implemented virtual base class methods virtual GClassCode code(void) const; - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual double mc_norm(const GSkyDir& dir, const double& radius) const; virtual void read(const GXmlElement& xml); virtual void write(GXmlElement& xml) const; diff --git a/include/GModelSpatialEllipticalDisk.hpp b/include/GModelSpatialEllipticalDisk.hpp index 77fad36..13717d0 100644 --- a/include/GModelSpatialEllipticalDisk.hpp +++ b/include/GModelSpatialEllipticalDisk.hpp @@ -68,11 +68,8 @@ public: virtual double eval(const double& theta, const double& posangle, const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const double& posangle, - const GEnergy& energy, - const GTime& time) const; + const GTime& time, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/include/GModelSpatialEllipticalGauss.hpp b/include/GModelSpatialEllipticalGauss.hpp index a2bc95f..3317ebf 100644 --- a/include/GModelSpatialEllipticalGauss.hpp +++ b/include/GModelSpatialEllipticalGauss.hpp @@ -69,11 +69,8 @@ public: virtual double eval(const double& theta, const double& posangle, const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const double& posangle, - const GEnergy& energy, - const GTime& time) const; + const GTime& time, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/include/GModelSpatialPointSource.hpp b/include/GModelSpatialPointSource.hpp index ba9049a..b5124e2 100644 --- a/include/GModelSpatialPointSource.hpp +++ b/include/GModelSpatialPointSource.hpp @@ -65,8 +65,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual GClassCode code(void) const; - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/include/GModelSpatialRadial.hpp b/include/GModelSpatialRadial.hpp index 3870c51..242af16 100644 --- a/include/GModelSpatialRadial.hpp +++ b/include/GModelSpatialRadial.hpp @@ -68,13 +68,11 @@ public: virtual std::string type(void) const = 0; virtual double eval(const double& theta, const GEnergy& energy, - const GTime& time) const = 0; - virtual double eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const = 0; + const GTime& time, + const bool& gradients = false) const = 0; virtual GSkyDir mc(const GEnergy& energy, - const GTime& time, - GRan& ran) const = 0; + const GTime& time, + GRan& ran) const = 0; virtual bool contains(const GSkyDir& dir, const double& margin = 0.0) const = 0; virtual double theta_max(void) const = 0; @@ -82,8 +80,8 @@ public: // Implemented pure virtual base class methods virtual GClassCode code(void) const; - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual double mc_norm(const GSkyDir& dir, const double& radius) const; virtual void read(const GXmlElement& xml); virtual void write(GXmlElement& xml) const; diff --git a/include/GModelSpatialRadialDisk.hpp b/include/GModelSpatialRadialDisk.hpp index 046aff3..d9bce74 100644 --- a/include/GModelSpatialRadialDisk.hpp +++ b/include/GModelSpatialRadialDisk.hpp @@ -65,10 +65,8 @@ public: virtual std::string type(void) const; virtual double eval(const double& theta, const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const; + const GTime& time, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/include/GModelSpatialRadialGauss.hpp b/include/GModelSpatialRadialGauss.hpp index 3ca8a3b..8bf6a17 100644 --- a/include/GModelSpatialRadialGauss.hpp +++ b/include/GModelSpatialRadialGauss.hpp @@ -66,10 +66,8 @@ public: virtual std::string type(void) const; virtual double eval(const double& theta, const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const; + const GTime& time, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/include/GModelSpatialRadialShell.hpp b/include/GModelSpatialRadialShell.hpp index 307fdf4..512c7ed 100644 --- a/include/GModelSpatialRadialShell.hpp +++ b/include/GModelSpatialRadialShell.hpp @@ -77,10 +77,8 @@ public: virtual std::string type(void) const; virtual double eval(const double& theta, const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const; + const GTime& time, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/include/GModelSpectral.hpp b/include/GModelSpectral.hpp index 32deb1e..3f3ca05 100644 --- a/include/GModelSpectral.hpp +++ b/include/GModelSpectral.hpp @@ -86,9 +86,8 @@ public: virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const = 0; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()) = 0; + const GTime& srcTime = GTime(), + const bool& gradients = false) const = 0; virtual double flux(const GEnergy& emin, const GEnergy& emax) const = 0; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralBrokenPlaw.hpp b/include/GModelSpectralBrokenPlaw.hpp index 25671cc..6a8e535 100644 --- a/include/GModelSpectralBrokenPlaw.hpp +++ b/include/GModelSpectralBrokenPlaw.hpp @@ -84,9 +84,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralConst.hpp b/include/GModelSpectralConst.hpp index 8faf21f..6f019f9 100644 --- a/include/GModelSpectralConst.hpp +++ b/include/GModelSpectralConst.hpp @@ -74,9 +74,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralExpInvPlaw.hpp b/include/GModelSpectralExpInvPlaw.hpp index fd7797c..20fa151 100644 --- a/include/GModelSpectralExpInvPlaw.hpp +++ b/include/GModelSpectralExpInvPlaw.hpp @@ -91,9 +91,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralExpPlaw.hpp b/include/GModelSpectralExpPlaw.hpp index 678e2c1..03bd4b8 100644 --- a/include/GModelSpectralExpPlaw.hpp +++ b/include/GModelSpectralExpPlaw.hpp @@ -87,9 +87,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralFunc.hpp b/include/GModelSpectralFunc.hpp index a30de1d..5c157f8 100644 --- a/include/GModelSpectralFunc.hpp +++ b/include/GModelSpectralFunc.hpp @@ -77,9 +77,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralGauss.hpp b/include/GModelSpectralGauss.hpp index 7a315ff..137b7ec 100644 --- a/include/GModelSpectralGauss.hpp +++ b/include/GModelSpectralGauss.hpp @@ -79,9 +79,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralLogParabola.hpp b/include/GModelSpectralLogParabola.hpp index 9a632b8..852c3be 100644 --- a/include/GModelSpectralLogParabola.hpp +++ b/include/GModelSpectralLogParabola.hpp @@ -88,9 +88,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralNodes.hpp b/include/GModelSpectralNodes.hpp index 5e2b0d8..78702b4 100644 --- a/include/GModelSpectralNodes.hpp +++ b/include/GModelSpectralNodes.hpp @@ -99,9 +99,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralPlaw.hpp b/include/GModelSpectralPlaw.hpp index ae81009..7373bd5 100644 --- a/include/GModelSpectralPlaw.hpp +++ b/include/GModelSpectralPlaw.hpp @@ -82,9 +82,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralPlawEnergyFlux.hpp b/include/GModelSpectralPlawEnergyFlux.hpp index 61103cd..32b759c 100644 --- a/include/GModelSpectralPlawEnergyFlux.hpp +++ b/include/GModelSpectralPlawEnergyFlux.hpp @@ -97,9 +97,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralPlawPhotonFlux.hpp b/include/GModelSpectralPlawPhotonFlux.hpp index f9ab16e..5b3b8e6 100644 --- a/include/GModelSpectralPlawPhotonFlux.hpp +++ b/include/GModelSpectralPlawPhotonFlux.hpp @@ -97,9 +97,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelSpectralSuperExpPlaw.hpp b/include/GModelSpectralSuperExpPlaw.hpp index cf8f5b0..9f3c701 100644 --- a/include/GModelSpectralSuperExpPlaw.hpp +++ b/include/GModelSpectralSuperExpPlaw.hpp @@ -92,9 +92,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/include/GModelTemporal.hpp b/include/GModelTemporal.hpp index efc5dba..48214c8 100644 --- a/include/GModelTemporal.hpp +++ b/include/GModelTemporal.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelTemporal.hpp - Abstract temporal model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2015 by Juergen Knoedlseder * + * copyright (C) 2009-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -68,8 +68,8 @@ public: virtual GModelTemporal* clone(void) const = 0; virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; - virtual double eval(const GTime& srcTime) const = 0; - virtual double eval_gradients(const GTime& srcTime) = 0; + virtual double eval(const GTime& srcTime, + const bool& gradients = false) const = 0; virtual GTimes mc(const double& rate, const GTime& tmin, const GTime& tmax, GRan& ran) const = 0; virtual void read(const GXmlElement& xml) = 0; diff --git a/include/GModelTemporalConst.hpp b/include/GModelTemporalConst.hpp index 2b42d01..f2bde66 100644 --- a/include/GModelTemporalConst.hpp +++ b/include/GModelTemporalConst.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelTemporalConst.hpp - Temporal constant model class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2014 by Juergen Knoedlseder * + * copyright (C) 2009-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -70,8 +70,8 @@ public: virtual GModelTemporalConst* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const GTime& srcTime) const; - virtual double eval_gradients(const GTime& srcTime); + virtual double eval(const GTime& srcTime, + const bool& gradients = false) const; virtual GTimes mc(const double& rate, const GTime& tmin, const GTime& tmax, GRan& ran) const; virtual void read(const GXmlElement& xml); diff --git a/include/GModels.hpp b/include/GModels.hpp index e42b774..4cd1c22 100644 --- a/include/GModels.hpp +++ b/include/GModels.hpp @@ -132,13 +132,17 @@ class GFilename; * models.read(xml); // Read models from GXml object * * The sum of all models in the container are evaluated for a given @p event - * and @p observation using the eval() or eval_gradients() methods: + * and @p observation using the eval() method: * * double value = models.eval(event, observation); - * double value = models.eval_gradients(event, observation); * - * The eval_gradients() method sets the parameter gradients for all free - * model parameters that have an analytical parameter gradient. + * If the eval() method is called with the optional gradients parameter set + * to true, i.e. + * + * double value = models.eval(event, observation, true); + * + * the method computes also the parameter gradients for all free model + * parameters that have an analytical parameter gradient. * * The only member of GModels is a list of model pointers. The class handles * the proper allocation and deallocation of the model memory. @@ -183,8 +187,9 @@ public: void write(GXml& xml) const; int npars(void) const; GOptimizerPars pars(void) const; - double eval(const GEvent& event, const GObservation& obs) const; - double eval_gradients(const GEvent& event, const GObservation& obs) const; + double eval(const GEvent& event, + const GObservation& obs, + const bool& gradients = false) const; std::string print(const GChatter& chatter = NORMAL) const; protected: diff --git a/include/GOptimizerPar.hpp b/include/GOptimizerPar.hpp index 587e700..e701ecc 100644 --- a/include/GOptimizerPar.hpp +++ b/include/GOptimizerPar.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GOptimizerPar.hpp - Optimizer parameter class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2015 by Juergen Knoedlseder * + * copyright (C) 2013-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -122,7 +122,7 @@ public: const double& scale(void) const; void factor_value(const double& value); void factor_error(const double& error); - void factor_gradient(const double& gradient); + void factor_gradient(const double& gradient) const; void factor_min(const double& min); void factor_max(const double& max); void factor_range(const double& min, const double& max); @@ -162,18 +162,18 @@ protected: void free_members(void); // Proteced data members - std::string m_name; //!< Parameter name - std::string m_unit; //!< Parameter unit - double m_factor_value; //!< Parameter value factor - double m_factor_error; //!< Uncertainty in parameter value factor - double m_factor_gradient; //!< Function gradient factor - double m_factor_min; //!< Parameter minimum factor - double m_factor_max; //!< Parameter maximum factor - double m_scale; //!< Parameter scaling (true = factor * scale) - bool m_free; //!< Parameter is free - bool m_has_min; //!< Parameter has minimum boundary - bool m_has_max; //!< Parameter has maximum boundary - bool m_has_grad; //!< Parameter has analytic gradient + std::string m_name; //!< Parameter name + std::string m_unit; //!< Parameter unit + double m_factor_value; //!< Parameter value factor + double m_factor_error; //!< Uncertainty in parameter value factor + double m_factor_min; //!< Parameter minimum factor + double m_factor_max; //!< Parameter maximum factor + mutable double m_factor_gradient; //!< Function gradient factor + double m_scale; //!< Parameter scaling (true = factor * scale) + bool m_free; //!< Parameter is free + bool m_has_min; //!< Parameter has minimum boundary + bool m_has_max; //!< Parameter has maximum boundary + bool m_has_grad; //!< Parameter has analytic gradient }; @@ -376,7 +376,7 @@ void GOptimizerPar::factor_error(const double& error) * Sets the parameter gradient factor. ***************************************************************************/ inline -void GOptimizerPar::factor_gradient(const double& gradient) +void GOptimizerPar::factor_gradient(const double& gradient) const { m_factor_gradient = gradient; return; diff --git a/inst/com/include/GCOMModelDRBFitting.hpp b/inst/com/include/GCOMModelDRBFitting.hpp index 66d39c9..e9e524a 100644 --- a/inst/com/include/GCOMModelDRBFitting.hpp +++ b/inst/com/include/GCOMModelDRBFitting.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCOMModelDRBFitting.hpp - COMPTEL DRB model fitting class * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2014 by Juergen Knoedlseder * + * copyright (C) 2012-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -67,9 +67,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; virtual GCOMEventCube* mc(const GObservation& obs, GRan& ran) const; diff --git a/inst/com/pyext/GCOMModelDRBFitting.i b/inst/com/pyext/GCOMModelDRBFitting.i index f404c56..9dd6a41 100644 --- a/inst/com/pyext/GCOMModelDRBFitting.i +++ b/inst/com/pyext/GCOMModelDRBFitting.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCOMModelDRBFitting.i - COMPTEL DRB model fitting class * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2015 by Juergen Knoedlseder * + * copyright (C) 2012-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -49,9 +49,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; virtual GCOMEventCube* mc(const GObservation& obs, GRan& ran) const; diff --git a/inst/com/src/GCOMModelDRBFitting.cpp b/inst/com/src/GCOMModelDRBFitting.cpp index 56761c9..5fa8102 100644 --- a/inst/com/src/GCOMModelDRBFitting.cpp +++ b/inst/com/src/GCOMModelDRBFitting.cpp @@ -43,9 +43,7 @@ const GCOMModelDRBFitting g_com_drb_fitting_seed; const GModelRegistry g_com_drb_fitting_registry(&g_com_drb_fitting_seed); /* __ Method name definitions ____________________________________________ */ -#define G_EVAL "GCOMModelDRBFitting::eval(GEvent&,GObservation&)" -#define G_EVAL_GRADIENTS "GCOMModelDRBFitting::eval_gradients(GEvent&,"\ - "GObservation&)" +#define G_EVAL "GCOMModelDRBFitting::eval(GEvent&, GObservation&, bool&)" #define G_READ "GCOMModelDRBFitting::read(GXmlElement&)" #define G_WRITE "GCOMModelDRBFitting::write(GXmlElement&)" @@ -229,6 +227,7 @@ GCOMModelDRBFitting* GCOMModelDRBFitting::clone(void) const * * @param[in] event Observed event. * @param[in] obs Observation. + * @param[in] gradients Compute gradients? * @return Background model value. * * @exception GException::invalid_argument @@ -238,7 +237,8 @@ GCOMModelDRBFitting* GCOMModelDRBFitting::clone(void) const * Evaluates the COMPTEL DRB fitting model. ***************************************************************************/ double GCOMModelDRBFitting::eval(const GEvent& event, - const GObservation& obs) const + const GObservation& obs, + const bool& gradients) const { // Extract COMPTEL observation const GCOMObservation* observation = dynamic_cast(&obs); @@ -262,107 +262,11 @@ double GCOMModelDRBFitting::eval(const GEvent& event, // Initialise value double value = 0.0; - // Get bin index - int index = bin->index(); - - // Get bin size. - double size = bin->size(); - - // Continue only if bin size is positive - if (size > 0.0) { - - // Initialise scaling factor - double scale = 1.0; - - // Get DRB model value - value = observation->drb().pixels()[index] / size; - - // If model is a scaling factor then use the single parameter as such - if (m_scale) { - scale = m_values[0].value(); - } - - // ... otherwise perform a linear interpolation - else { - - // Get Phibar value - double phibar = bin->dir().phibar(); - - // Update evaluation cache - update_cache(); - - // Set node array for linear interpolation - m_nodes.set_value(phibar); - - // Get scale factor - scale = m_values[m_nodes.inx_left()].value() * m_nodes.wgt_left() + - m_values[m_nodes.inx_right()].value() * m_nodes.wgt_right(); - - } // endelse: performed linear interpolation - - // Compute background value - value *= scale; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GCOMModelDRBFitting::eval"; - std::cout << "(index=" << index << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", scale=" << scale; - std::cout << ")" << std::endl; + // Optionally initialise gradients + if (gradients) { + for (int i = 0; i < m_values.size(); ++i) { + m_values[i].factor_gradient(0.0); } - #endif - - } // endif: binsize was positive - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] event Observed event. - * @param[in] obs Observation. - * @return Background model value. - * - * @exception GCOMException::bad_observation_type - * Observation is not a COMPTEL observation. - * @exception GCOMException::bad_event_type - * Event is not a COMPTEL event bin. - * - * Evaluates the COMPTEL DRB fitting model value and sets the parameter - * gradients. - ***************************************************************************/ -double GCOMModelDRBFitting::eval_gradients(const GEvent& event, - const GObservation& obs) const -{ - // Extract COMPTEL observation - const GCOMObservation* observation = dynamic_cast(&obs); - if (observation == NULL) { - std::string cls = std::string(typeid(&obs).name()); - std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL " - "observations. Please specify a COMPTEL observation " - "as argument."; - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Extract COMPTEL event bin - const GCOMEventBin* bin = dynamic_cast(&event); - if (bin == NULL) { - std::string cls = std::string(typeid(&event).name()); - std::string msg = "Event of type \""+cls+"\" is not a COMPTEL event. " - "Please specify a COMPTEL event as argument."; - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Initialise value and gradients - double value = 0.0; - for (int i = 0; i < m_values.size(); ++i) { - const_cast(this)->m_values[i].factor_gradient(0.0); } // Get bin index @@ -380,20 +284,26 @@ double GCOMModelDRBFitting::eval_gradients(const GEvent& event, // Get DRB model value value = observation->drb().pixels()[index] / size; - - // If model is a scaling factor then use the single parameter as such + // If the model has a single scaling factor then use the first + // parameter as scaling factor of the model if (m_scale) { // Get scale factor scale = m_values[0].value(); - - // Compute partial derivative - double grad = (m_values[0].is_free()) ? value * m_values[0].scale() : 0.0; - // Set gradient (circumvent const correctness) - const_cast(this)->m_values[0].factor_gradient(grad); + // Optionally compute gradients + if (gradients) { - } + // Compute partial derivative + double grad = (m_values[0].is_free()) + ? value * m_values[0].scale() : 0.0; + + // Set gradient + m_values[0].factor_gradient(grad); + + } // endif: gradient computation was requested + + } // endif: model is a scaling factor // ... otherwise perform a linear interpolation else { @@ -417,17 +327,22 @@ double GCOMModelDRBFitting::eval_gradients(const GEvent& event, scale = m_values[inx_left].value() * wgt_left + m_values[inx_right].value() * wgt_right; - // Gradient for left node - if (m_values[inx_left].is_free()) { - double grad = wgt_left * value * m_values[inx_left].scale(); - const_cast(this)->m_values[inx_left].factor_gradient(grad); - } + // Optionally compute gradients + if (gradients) { - // Gradient for right node - if (m_values[inx_right].is_free()) { - double grad = wgt_right * value * m_values[inx_right].scale(); - const_cast(this)->m_values[inx_right].factor_gradient(grad); - } + // Compute partial derivatives + double g_left = (m_values[inx_left].is_free()) + ? wgt_left * value * m_values[inx_left].scale() + : 0.0; + double g_right = (m_values[inx_right].is_free()) + ? wgt_right * value * m_values[inx_right].scale() + : 0.0; + + // Set gradients + m_values[inx_left].factor_gradient(g_left); + m_values[inx_right].factor_gradient(g_right); + + } // endif: gradient computation was requested } // endelse: performed linear interpolation @@ -437,7 +352,7 @@ double GCOMModelDRBFitting::eval_gradients(const GEvent& event, // Compile option: Check for NaN/Inf #if defined(G_NAN_CHECK) if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GCOMModelDRBFitting::eval_gradients"; + std::cout << "*** ERROR: GCOMModelDRBFitting::eval"; std::cout << "(index=" << index << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; diff --git a/inst/cta/include/GCTAModelAeffBackground.hpp b/inst/cta/include/GCTAModelAeffBackground.hpp index ed65c27..05ca443 100644 --- a/inst/cta/include/GCTAModelAeffBackground.hpp +++ b/inst/cta/include/GCTAModelAeffBackground.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelAeffBackground.hpp - CTA Aeff background model class * * ----------------------------------------------------------------------- * - * copyright (C) 2015 by Michael Mayer * + * copyright (C) 2015-2016 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -64,9 +64,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; diff --git a/inst/cta/include/GCTAModelCubeBackground.hpp b/inst/cta/include/GCTAModelCubeBackground.hpp index 25f05bf..7da5052 100644 --- a/inst/cta/include/GCTAModelCubeBackground.hpp +++ b/inst/cta/include/GCTAModelCubeBackground.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelCubeBackground.hpp - CTA cube background model class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2015 by Michael Mayer * + * copyright (C) 2013-2016 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -69,9 +69,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; diff --git a/inst/cta/include/GCTAModelIrfBackground.hpp b/inst/cta/include/GCTAModelIrfBackground.hpp index ddb56fc..08e86a1 100644 --- a/inst/cta/include/GCTAModelIrfBackground.hpp +++ b/inst/cta/include/GCTAModelIrfBackground.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelIrfBackground.hpp - CTA IRF background model class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2015 by Juergen Knoedlseder * + * copyright (C) 2014-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -64,9 +64,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; diff --git a/inst/cta/include/GCTAModelRadial.hpp b/inst/cta/include/GCTAModelRadial.hpp index cf29afe..f881b28 100644 --- a/inst/cta/include/GCTAModelRadial.hpp +++ b/inst/cta/include/GCTAModelRadial.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadial.hpp - Radial model abstract base class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -64,8 +64,8 @@ public: virtual GCTAModelRadial* clone(void) const = 0; virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; - virtual double eval(const double& offset) const = 0; - virtual double eval_gradients(const double& offset) const = 0; + virtual double eval(const double& offset, + const bool& gradients = false) const = 0; virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const = 0; virtual double omega(void) const = 0; virtual void read(const GXmlElement& xml) = 0; diff --git a/inst/cta/include/GCTAModelRadialAcceptance.hpp b/inst/cta/include/GCTAModelRadialAcceptance.hpp index d94095e..a231274 100644 --- a/inst/cta/include/GCTAModelRadialAcceptance.hpp +++ b/inst/cta/include/GCTAModelRadialAcceptance.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialAcceptance.hpp - Radial acceptance model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2014 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -70,9 +70,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const; diff --git a/inst/cta/include/GCTAModelRadialGauss.hpp b/inst/cta/include/GCTAModelRadialGauss.hpp index 74ec039..b0b37d0 100644 --- a/inst/cta/include/GCTAModelRadialGauss.hpp +++ b/inst/cta/include/GCTAModelRadialGauss.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialGauss.hpp - Radial Gaussian CTA model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2014 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -69,8 +69,8 @@ public: virtual GCTAModelRadialGauss* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const double& offset) const; - virtual double eval_gradients(const double& offset) const; + virtual double eval(const double& offset, + const bool& gradients = false) const; virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const; virtual double omega(void) const; virtual void read(const GXmlElement& xml); @@ -78,8 +78,8 @@ public: virtual std::string print(const GChatter& chatter = NORMAL) const; // Other methods - double sigma(void) const { return m_sigma.value(); } - void sigma(const double& sigma) { m_sigma.value(sigma); } + double sigma(void) const; + void sigma(const double& sigma); protected: // Protected methods @@ -130,4 +130,29 @@ std::string GCTAModelRadialGauss::type(void) const return ("Gaussian"); } + +/***********************************************************************//** + * @brief Return Gaussian width parameter + * + * @return Gaussian width parameter. + ***************************************************************************/ +inline +double GCTAModelRadialGauss::sigma(void) const +{ + return (m_sigma.value()); +} + + +/***********************************************************************//** + * @brief Set Gaussian width parameter + * + * @param[in] sigma Gaussian width parameter. + ***************************************************************************/ +inline +void GCTAModelRadialGauss::sigma(const double& sigma) +{ + m_sigma.value(sigma); + return; +} + #endif /* GCTAMODELRADIALGAUSS_HPP */ diff --git a/inst/cta/include/GCTAModelRadialPolynom.hpp b/inst/cta/include/GCTAModelRadialPolynom.hpp index ebfcf8b..b3111f4 100644 --- a/inst/cta/include/GCTAModelRadialPolynom.hpp +++ b/inst/cta/include/GCTAModelRadialPolynom.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialPolynom.hpp - Radial Polynom CTA model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -70,8 +70,8 @@ public: virtual GCTAModelRadialPolynom* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const double& offset) const; - virtual double eval_gradients(const double& offset) const; + virtual double eval(const double& offset, + const bool& gradients = false) const; virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const; virtual double omega(void) const; virtual void read(const GXmlElement& xml); diff --git a/inst/cta/include/GCTAModelRadialProfile.hpp b/inst/cta/include/GCTAModelRadialProfile.hpp index e2729b1..95e678a 100644 --- a/inst/cta/include/GCTAModelRadialProfile.hpp +++ b/inst/cta/include/GCTAModelRadialProfile.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialProfile.hpp - Radial Profile CTA model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2014 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -71,8 +71,8 @@ public: virtual GCTAModelRadialProfile* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const double& offset) const; - virtual double eval_gradients(const double& offset) const; + virtual double eval(const double& offset, + const bool& gradients = false) const; virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const; virtual double omega(void) const; virtual void read(const GXmlElement& xml); diff --git a/inst/cta/pyext/GCTAModelAeffBackground.i b/inst/cta/pyext/GCTAModelAeffBackground.i index f222238..13b5fbe 100644 --- a/inst/cta/pyext/GCTAModelAeffBackground.i +++ b/inst/cta/pyext/GCTAModelAeffBackground.i @@ -1,7 +1,7 @@ /*************************************************************************** - * GCTAModelAeffBackground.i - CTA IRF background model class * + * GCTAModelAeffBackground.i - CTA Aeff background model class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2015 by Juergen Knoedlseder * + * copyright (C) 2015-2016 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -44,21 +44,20 @@ public: virtual ~GCTAModelAeffBackground(void); // Implemented pure virtual methods - virtual void clear(void); + virtual void clear(void); virtual GCTAModelAeffBackground* clone(void) const; - virtual std::string classname(void) const; - virtual std::string type(void) const; - virtual bool is_constant(void) const; - virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; - virtual double npred(const GEnergy& obsEng, - const GTime& obsTime, - const GObservation& obs) const; - virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const; - virtual void read(const GXmlElement& xml); - virtual void write(GXmlElement& xml) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual bool is_constant(void) const; + virtual double eval(const GEvent& event, + const GObservation& obs, + const bool& gradients = false) const; + virtual double npred(const GEnergy& obsEng, + const GTime& obsTime, + const GObservation& obs) const; + virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; // Other methods GModelSpectral* spectral(void) const; diff --git a/inst/cta/pyext/GCTAModelCubeBackground.i b/inst/cta/pyext/GCTAModelCubeBackground.i index dccaea9..ce1cdf1 100644 --- a/inst/cta/pyext/GCTAModelCubeBackground.i +++ b/inst/cta/pyext/GCTAModelCubeBackground.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelCubeBackground.i - CTA cube background model class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2015 by Michael Mayer * + * copyright (C) 2013-2016 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -52,9 +52,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const; diff --git a/inst/cta/pyext/GCTAModelIrfBackground.i b/inst/cta/pyext/GCTAModelIrfBackground.i index f50b9da..a75f4d6 100644 --- a/inst/cta/pyext/GCTAModelIrfBackground.i +++ b/inst/cta/pyext/GCTAModelIrfBackground.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelIrfBackground.i - CTA IRF background model class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2015 by Juergen Knoedlseder * + * copyright (C) 2014-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -50,9 +50,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; diff --git a/inst/cta/pyext/GCTAModelRadial.i b/inst/cta/pyext/GCTAModelRadial.i index 7997645..ec8a7b7 100644 --- a/inst/cta/pyext/GCTAModelRadial.i +++ b/inst/cta/pyext/GCTAModelRadial.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadial.i - Abstract radial model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -48,8 +48,8 @@ public: virtual GCTAModelRadial* clone(void) const = 0; virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; - virtual double eval(const double& offset) const = 0; - virtual double eval_gradients(const double& offset) const = 0; + virtual double eval(const double& offset, + const bool& gradients = false) const = 0; virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const = 0; virtual double omega(void) const = 0; virtual void read(const GXmlElement& xml) = 0; diff --git a/inst/cta/pyext/GCTAModelRadialAcceptance.i b/inst/cta/pyext/GCTAModelRadialAcceptance.i index fe80308..3d4dc65 100644 --- a/inst/cta/pyext/GCTAModelRadialAcceptance.i +++ b/inst/cta/pyext/GCTAModelRadialAcceptance.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialAcceptance.i - Radial acceptance model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -53,9 +53,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const; diff --git a/inst/cta/pyext/GCTAModelRadialGauss.i b/inst/cta/pyext/GCTAModelRadialGauss.i index e195584..8f195db 100644 --- a/inst/cta/pyext/GCTAModelRadialGauss.i +++ b/inst/cta/pyext/GCTAModelRadialGauss.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialGauss.i - Radial Gaussian model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -49,8 +49,8 @@ public: virtual GCTAModelRadialGauss* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const double& offset) const; - virtual double eval_gradients(const double& offset) const; + virtual double eval(const double& offset, + const bool& gradients = false) const; virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const; virtual double omega(void) const; virtual void read(const GXmlElement& xml); diff --git a/inst/cta/pyext/GCTAModelRadialPolynom.i b/inst/cta/pyext/GCTAModelRadialPolynom.i index dc983d5..90d89bd 100644 --- a/inst/cta/pyext/GCTAModelRadialPolynom.i +++ b/inst/cta/pyext/GCTAModelRadialPolynom.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialPolynom.i - Radial Polynom model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -49,8 +49,8 @@ public: virtual GCTAModelRadialPolynom* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const double& offset) const; - virtual double eval_gradients(const double& offset) const; + virtual double eval(const double& offset, + const bool& gradients = false) const; virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const; virtual double omega(void) const; virtual void read(const GXmlElement& xml); diff --git a/inst/cta/pyext/GCTAModelRadialProfile.i b/inst/cta/pyext/GCTAModelRadialProfile.i index 1f8c8d9..44f9dba 100644 --- a/inst/cta/pyext/GCTAModelRadialProfile.i +++ b/inst/cta/pyext/GCTAModelRadialProfile.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialProfile.i - Radial Profile model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -50,8 +50,8 @@ public: virtual GCTAModelRadialProfile* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const double& offset) const; - virtual double eval_gradients(const double& offset) const; + virtual double eval(const double& offset, + const bool& gradients = false) const; virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const; virtual double omega(void) const; virtual void read(const GXmlElement& xml); diff --git a/inst/cta/src/GCTAModelAeffBackground.cpp b/inst/cta/src/GCTAModelAeffBackground.cpp index ccb4daa..05594c2 100644 --- a/inst/cta/src/GCTAModelAeffBackground.cpp +++ b/inst/cta/src/GCTAModelAeffBackground.cpp @@ -46,9 +46,7 @@ const GCTAModelAeffBackground g_cta_aeff_background_seed; const GModelRegistry g_cta_aeff_background_registry(&g_cta_aeff_background_seed); /* __ Method name definitions ____________________________________________ */ -#define G_EVAL "GCTAModelAeffBackground::eval(GEvent&, GObservation&)" -#define G_EVAL_GRADIENTS "GCTAModelAeffBackground::eval_gradients(GEvent&,"\ - " GObservation&)" +#define G_EVAL "GCTAModelAeffBackground::eval(GEvent&, GObservation&, bool&)" #define G_NPRED "GCTAModelAeffBackground::npred(GEnergy&, GTime&,"\ " GObservation&)" #define G_MC "GCTAModelAeffBackground::mc(GObservation&, GRan&)" @@ -274,15 +272,20 @@ GCTAModelAeffBackground* GCTAModelAeffBackground::clone(void) const * * @param[in] event Observed event. * @param[in] obs Observation. + * @param[in] gradients Compute gradients? * @return Function value. * * @exception GException::invalid_argument * Specified observation is not of the expected type. * + * If the @p gradients flag is true the method will also set the parameter + * gradients of the model parameters. + * * @todo Make sure that DETX and DETY are always set in GCTAInstDir. ***************************************************************************/ double GCTAModelAeffBackground::eval(const GEvent& event, - const GObservation& obs) const + const GObservation& obs, + const bool& gradients) const { // Get pointer on CTA observation const GCTAObservation* cta = dynamic_cast(&obs); @@ -330,85 +333,10 @@ double GCTAModelAeffBackground::eval(const GEvent& event, cta->pointing().zenith(), cta->pointing().azimuth(), false); double spec = (spectral() != NULL) - ? spectral()->eval(event.energy(), event.time()) : 1.0; + ? spectral()->eval(event.energy(), event.time(), gradients) + : 1.0; double temp = (temporal() != NULL) - ? temporal()->eval(event.time()) : 1.0; - - // Compute value - double value = spat * spec * temp; - - // Apply deadtime correction - value *= obs.deadc(event.time()); - - // Return value - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] event Observed event. - * @param[in] obs Observation. - * @return Function value. - * - * @exception GException::invalid_argument - * Specified observation is not of the expected type. - * - * @todo Make sure that DETX and DETY are always set in GCTAInstDir. - ***************************************************************************/ -double GCTAModelAeffBackground::eval_gradients(const GEvent& event, - const GObservation& obs) const -{ - // Get pointer on CTA observation - const GCTAObservation* cta = dynamic_cast(&obs); - if (cta == NULL) { - std::string msg = "Specified observation is not a CTA observation.\n" + - obs.print(); - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Get pointer on CTA IRF response - const GCTAResponseIrf* rsp = dynamic_cast(cta->response()); - if (rsp == NULL) { - std::string msg = "Specified observation does not contain an" - " IRF response.\n" + obs.print(); - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Retrieve pointer to CTA Effective Area - const GCTAAeff* aeff = rsp->aeff(); - if (aeff == NULL) { - std::string msg = "Specified observation contains no effective area" - " information.\n" + obs.print(); - throw GException::invalid_argument(G_EVAL, msg); - } - - // Extract CTA instrument direction from event - const GCTAInstDir* dir = dynamic_cast(&(event.dir())); - if (dir == NULL) { - std::string msg = "No CTA instrument direction found in event."; - throw GException::invalid_argument(G_EVAL, msg); - } - - // Set DETX and DETY in instrument direction - GCTAInstDir inst_dir = cta->pointing().instdir(dir->dir()); - - // Set theta and phi from instrument coordinates - double theta = std::sqrt(inst_dir.detx() * inst_dir.detx() + - inst_dir.dety() * inst_dir.dety()); - double phi = gammalib::atan2d(inst_dir.dety(), inst_dir.detx()) * - gammalib::deg2rad; - - // Evaluate function - double logE = event.energy().log10TeV(); - double spat = (*aeff)(logE, theta, phi, - cta->pointing().zenith(), - cta->pointing().azimuth(), false); - double spec = (spectral() != NULL) - ? spectral()->eval_gradients(event.energy(), event.time()) : 1.0; - double temp = (temporal() != NULL) - ? temporal()->eval_gradients(event.time()) : 1.0; + ? temporal()->eval(event.time(), gradients) : 1.0; // Compute value double value = spat * spec * temp; @@ -417,23 +345,28 @@ double GCTAModelAeffBackground::eval_gradients(const GEvent& event, double deadc = obs.deadc(event.time()); value *= deadc; - // Multiply factors to spectral gradients - if (spectral() != NULL) { - double fact = spat * temp * deadc; - if (fact != 1.0) { - for (int i = 0; i < spectral()->size(); ++i) - (*spectral())[i].factor_gradient( (*spectral())[i].factor_gradient() * fact ); + // Optionally compute partial derivatives + if (gradients) { + + // Multiply factors to spectral gradients + if (spectral() != NULL) { + double fact = spat * temp * deadc; + if (fact != 1.0) { + for (int i = 0; i < spectral()->size(); ++i) + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact ); + } } - } - // Multiply factors to temporal gradients - if (temporal() != NULL) { - double fact = spat * spec * deadc; - if (fact != 1.0) { - for (int i = 0; i < temporal()->size(); ++i) - (*temporal())[i].factor_gradient( (*temporal())[i].factor_gradient() * fact ); + // Multiply factors to temporal gradients + if (temporal() != NULL) { + double fact = spat * spec * deadc; + if (fact != 1.0) { + for (int i = 0; i < temporal()->size(); ++i) + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact ); + } } - } + + } // endif: computed partial derivatives // Return value return value; diff --git a/inst/cta/src/GCTAModelCubeBackground.cpp b/inst/cta/src/GCTAModelCubeBackground.cpp index d8b1715..4051716 100644 --- a/inst/cta/src/GCTAModelCubeBackground.cpp +++ b/inst/cta/src/GCTAModelCubeBackground.cpp @@ -44,9 +44,7 @@ const GCTAModelCubeBackground g_cta_inst_background_seed; const GModelRegistry g_cta_inst_background_registry(&g_cta_inst_background_seed); /* __ Method name definitions ____________________________________________ */ -#define G_EVAL "GCTAModelCubeBackground::eval(GEvent&, GObservation&)" -#define G_EVAL_GRADIENTS "GCTAModelCubeBackground::eval_gradients(GEvent&,"\ - " GObservation&)" +#define G_EVAL "GCTAModelCubeBackground::eval(GEvent&, GObservation&, bool&)" #define G_NPRED "GCTAModelCubeBackground::npred(GEnergy&, GTime&,"\ " GObservation&)" #define G_MC "GCTAModelCubeBackground::mc(GObservation&, GRan&)" @@ -251,13 +249,18 @@ GCTAModelCubeBackground* GCTAModelCubeBackground::clone(void) const * * @param[in] event Observed event. * @param[in] obs Observation. + * @param[in] gradients Compute gradients? * @return Function value. * * @exception GException::invalid_argument * Specified observation is not of the expected type. + * + * If the @p gradients flag is true the method will also set the parameter + * gradients of the model parameters. ***************************************************************************/ double GCTAModelCubeBackground::eval(const GEvent& event, - const GObservation& obs) const + const GObservation& obs, + const bool& gradients) const { // Get pointer on CTA observation const GCTAObservation* cta = dynamic_cast(&obs); @@ -287,89 +290,38 @@ double GCTAModelCubeBackground::eval(const GEvent& event, //double logE = event.energy().log10TeV(); double spat = bgd((*dir), event.energy()); double spec = (spectral() != NULL) - ? spectral()->eval(event.energy(), event.time()) : 1.0; - double temp = (temporal() != NULL) - ? temporal()->eval(event.time()) : 1.0; - - // Compute value. Note that background rates are already per - // livetime, hence no deadtime correction is needed here. - double value = spat * spec * temp; - - // Return value - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] event Observed event. - * @param[in] obs Observation. - * @return Function value. - * - * @exception GException::invalid_argument - * Specified observation is not of the expected type. - ***************************************************************************/ -double GCTAModelCubeBackground::eval_gradients(const GEvent& event, - const GObservation& obs) const -{ - // Get pointer on CTA observation - const GCTAObservation* cta = dynamic_cast(&obs); - if (cta == NULL) { - std::string msg = "Specified observation is not a CTA observation.\n" + - obs.print(); - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Get pointer on CTA IRF response - const GCTAResponseCube* rsp = dynamic_cast(cta->response()); - if (rsp == NULL) { - std::string msg = "Specified observation does not contain a" - " cube response.\n" + obs.print(); - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Extract CTA instrument direction from event - const GCTAInstDir* dir = dynamic_cast(&(event.dir())); - if (dir == NULL) { - std::string msg = "No CTA instrument direction found in event."; - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Retrieve reference to CTA cube background - const GCTACubeBackground& bgd = rsp->background(); - - // Evaluate function - //double logE = event.energy().log10TeV(); - double spat = bgd((*dir), event.energy()); - double spec = (spectral() != NULL) - ? spectral()->eval_gradients(event.energy(), event.time()) + ? spectral()->eval(event.energy(), event.time(), gradients) : 1.0; double temp = (temporal() != NULL) - ? temporal()->eval_gradients(event.time()) - : 1.0; + ? temporal()->eval(event.time(), gradients) : 1.0; // Compute value. Note that background rates are already per // livetime, hence no deadtime correction is needed here. double value = spat * spec * temp; - // Multiply factors to spectral gradients - if (spectral() != NULL) { - double fact = spat * temp; - if (fact != 1.0) { - for (int i = 0; i < spectral()->size(); ++i) - (*spectral())[i].factor_gradient( (*spectral())[i].factor_gradient() * fact ); + // If gradients were requested then multiply factors to spectral and + // temporal gradients + if (gradients) { + + // Multiply factors to spectral gradients + if (spectral() != NULL) { + double fact = spat * temp; + if (fact != 1.0) { + for (int i = 0; i < spectral()->size(); ++i) + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact ); + } } - } - // Multiply factors to temporal gradients - if (temporal() != NULL) { - double fact = spat * spec; - if (fact != 1.0) { - for (int i = 0; i < temporal()->size(); ++i) - (*temporal())[i].factor_gradient( (*temporal())[i].factor_gradient() * fact ); + // Multiply factors to temporal gradients + if (temporal() != NULL) { + double fact = spat * spec; + if (fact != 1.0) { + for (int i = 0; i < temporal()->size(); ++i) + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact ); + } } - } + + } // endif: gradients were requested // Return value return value; diff --git a/inst/cta/src/GCTAModelIrfBackground.cpp b/inst/cta/src/GCTAModelIrfBackground.cpp index 58c9ded..1e74735 100644 --- a/inst/cta/src/GCTAModelIrfBackground.cpp +++ b/inst/cta/src/GCTAModelIrfBackground.cpp @@ -46,9 +46,7 @@ const GCTAModelIrfBackground g_cta_inst_background_seed; const GModelRegistry g_cta_inst_background_registry(&g_cta_inst_background_seed); /* __ Method name definitions ____________________________________________ */ -#define G_EVAL "GCTAModelIrfBackground::eval(GEvent&, GObservation&)" -#define G_EVAL_GRADIENTS "GCTAModelIrfBackground::eval_gradients(GEvent&,"\ - " GObservation&)" +#define G_EVAL "GCTAModelIrfBackground::eval(GEvent&, GObservation&, bool&)" #define G_NPRED "GCTAModelIrfBackground::npred(GEnergy&, GTime&,"\ " GObservation&)" #define G_MC "GCTAModelIrfBackground::mc(GObservation&, GRan&)" @@ -274,15 +272,20 @@ GCTAModelIrfBackground* GCTAModelIrfBackground::clone(void) const * * @param[in] event Observed event. * @param[in] obs Observation. + * @param[in] gradients Compute gradients? * @return Function value. * * @exception GException::invalid_argument * Specified observation is not of the expected type. * + * If the @p gradients flag is true the method will also set the parameter + * gradients of the model parameters. + * * @todo Make sure that DETX and DETY are always set in GCTAInstDir. ***************************************************************************/ -double GCTAModelIrfBackground::eval(const GEvent& event, - const GObservation& obs) const +double GCTAModelIrfBackground::eval(const GEvent& event, + const GObservation& obs, + const bool& gradients) const { // Get pointer on CTA observation const GCTAObservation* cta = dynamic_cast(&obs); @@ -322,79 +325,10 @@ double GCTAModelIrfBackground::eval(const GEvent& event, double logE = event.energy().log10TeV(); double spat = (*bgd)(logE, inst_dir.detx(), inst_dir.dety()); double spec = (spectral() != NULL) - ? spectral()->eval(event.energy(), event.time()) : 1.0; - double temp = (temporal() != NULL) - ? temporal()->eval(event.time()) : 1.0; - - // Compute value - double value = spat * spec * temp; - - // Apply deadtime correction - value *= obs.deadc(event.time()); - - // Return value - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] event Observed event. - * @param[in] obs Observation. - * @return Function value. - * - * @exception GException::invalid_argument - * Specified observation is not of the expected type. - * - * @todo Make sure that DETX and DETY are always set in GCTAInstDir. - ***************************************************************************/ -double GCTAModelIrfBackground::eval_gradients(const GEvent& event, - const GObservation& obs) const -{ - // Get pointer on CTA observation - const GCTAObservation* cta = dynamic_cast(&obs); - if (cta == NULL) { - std::string msg = "Specified observation is not a CTA observation.\n" + - obs.print(); - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Get pointer on CTA IRF response - const GCTAResponseIrf* rsp = dynamic_cast(cta->response()); - if (rsp == NULL) { - std::string msg = "Specified observation does not contain an" - " IRF response.\n" + obs.print(); - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Retrieve pointer to CTA background - const GCTABackground* bgd = rsp->background(); - if (bgd == NULL) { - std::string msg = "Specified observation contains no background" - " information.\n" + obs.print(); - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Extract CTA instrument direction from event - const GCTAInstDir* dir = dynamic_cast(&(event.dir())); - if (dir == NULL) { - std::string msg = "No CTA instrument direction found in event."; - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); - } - - // Set DETX and DETY in instrument direction - GCTAInstDir inst_dir = cta->pointing().instdir(dir->dir()); - - // Evaluate function - double logE = event.energy().log10TeV(); - double spat = (*bgd)(logE, inst_dir.detx(), inst_dir.dety()); - double spec = (spectral() != NULL) - ? spectral()->eval_gradients(event.energy(), event.time()) + ? spectral()->eval(event.energy(), event.time(), gradients) : 1.0; double temp = (temporal() != NULL) - ? temporal()->eval_gradients(event.time()) - : 1.0; + ? temporal()->eval(event.time(), gradients) : 1.0; // Compute value double value = spat * spec * temp; @@ -403,23 +337,29 @@ double GCTAModelIrfBackground::eval_gradients(const GEvent& event, double deadc = obs.deadc(event.time()); value *= deadc; - // Multiply factors to spectral gradients - if (spectral() != NULL) { - double fact = spat * temp * deadc; - if (fact != 1.0) { - for (int i = 0; i < spectral()->size(); ++i) - (*spectral())[i].factor_gradient( (*spectral())[i].factor_gradient() * fact ); + // If gradients were requested then multiply factors to spectral and + // temporal gradients + if (gradients) { + + // Multiply factors to spectral gradients + if (spectral() != NULL) { + double fact = spat * temp * deadc; + if (fact != 1.0) { + for (int i = 0; i < spectral()->size(); ++i) + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact ); + } } - } - // Multiply factors to temporal gradients - if (temporal() != NULL) { - double fact = spat * spec * deadc; - if (fact != 1.0) { - for (int i = 0; i < temporal()->size(); ++i) - (*temporal())[i].factor_gradient( (*temporal())[i].factor_gradient() * fact ); + // Multiply factors to temporal gradients + if (temporal() != NULL) { + double fact = spat * spec * deadc; + if (fact != 1.0) { + for (int i = 0; i < temporal()->size(); ++i) + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact ); + } } - } + + } // endif: gradients were requested // Return value return value; diff --git a/inst/cta/src/GCTAModelRadialAcceptance.cpp b/inst/cta/src/GCTAModelRadialAcceptance.cpp index ecea7ba..268209a 100644 --- a/inst/cta/src/GCTAModelRadialAcceptance.cpp +++ b/inst/cta/src/GCTAModelRadialAcceptance.cpp @@ -53,10 +53,8 @@ const GModelRegistry g_cta_radial_acceptance_registry(&g_cta_radial_a /* __ Method name definitions ____________________________________________ */ #define G_ACCESS "GCTAModelRadialAcceptance::operator() (int)" -#define G_EVAL "GCTAModelRadialAcceptance::eval(GEvent&,"\ - " GObservation&)" -#define G_EVAL_GRADIENTS "GCTAModelRadialAcceptance::eval_gradients(GEvent&,"\ - " GObservation&)" +#define G_EVAL "GCTAModelRadialAcceptance::eval(GEvent&, GObservation&, "\ + "bool&)" #define G_NPRED "GCTAModelRadialAcceptance::npred(GEnergy&, GTime&,"\ " GObservation&)" #define G_MC "GCTAModelRadialAcceptance::mc(GObservation&, GRan&)" @@ -282,77 +280,18 @@ GCTAModelRadialAcceptance* GCTAModelRadialAcceptance::clone(void) const * * @param[in] event Observed event. * @param[in] obs Observation. + * @param[in] gradients Compute gradients? + * @return Function value. * * Evaluates tha CTA radial acceptance model which is a factorization of a * spatial, spectral and temporal model component. This method also applies * a deadtime correction factor, so that the normalization of the model is * a real rate (counts/exposure time). * - * @todo Add bookkeeping of last value and evaluate only if argument - * changed - * @todo Verify that CTA instrument direction pointer is valid, or better, - * add an offset method to GCTAPointing. Ideally, we should precompute - * all offset angles (for an event cube this may only be done easily - * if the pointing has been fixed; otherwise we need a structure - * similar to the Fermi/LAT livetime cube that provides the effective - * sky exposure as function of offset angle). - ***************************************************************************/ -double GCTAModelRadialAcceptance::eval(const GEvent& event, - const GObservation& obs) const -{ - // Get pointer on CTA observation - const GCTAObservation* ctaobs = dynamic_cast(&obs); - if (ctaobs == NULL) { - std::string msg = "Specified observation is not a CTA observation.\n" + - obs.print(); - throw GException::invalid_argument(G_EVAL, msg); - } - - // Get pointer on CTA pointing - const GCTAPointing& pnt = ctaobs->pointing(); - - // Get instrument direction - const GInstDir* inst_dir = &(event.dir()); - const GCTAInstDir* cta_dir = static_cast(inst_dir); - - // Compute offset angle (in degrees) - double offset = cta_dir->dir().dist_deg(pnt.dir()); - - // Evaluate function and gradients - double rad = (radial() != NULL) - ? radial()->eval(offset) : 1.0; - double spec = (spectral() != NULL) - ? spectral()->eval(event.energy(), event.time()) : 1.0; - double temp = (temporal() != NULL) - ? temporal()->eval(event.time()) : 1.0; - - // Compute value - double value = rad * spec * temp; - - // Apply deadtime correction - value *= obs.deadc(event.time()); - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] event Observed event. - * @param[in] obs Observation (not used). + * If the @p gradients flag is true the method will also set the parameter + * gradients of the model parameters. * - * @exception GCTAException::no_pointing - * No valid CTA pointing found in observation - * - * Evaluates tha CTA radial acceptance model and parameter gradients. The CTA - * radial acceptance model is a factorization of a spatial, spectral and - * temporal model component. This method also applies a deadtime correction - * factor, so that the normalization of the model is a real rate - * (counts/exposure time). - * - * @todo Add bookkeeping of last value and evaluate only if argument + * @todo Add bookkeeping of last value and evaluate only if argument * changed * @todo Verify that CTA instrument direction pointer is valid, or better, * add an offset method to GCTAPointing. Ideally, we should precompute @@ -361,15 +300,16 @@ double GCTAModelRadialAcceptance::eval(const GEvent& event, * similar to the Fermi/LAT livetime cube that provides the effective * sky exposure as function of offset angle). ***************************************************************************/ -double GCTAModelRadialAcceptance::eval_gradients(const GEvent& event, - const GObservation& obs) const +double GCTAModelRadialAcceptance::eval(const GEvent& event, + const GObservation& obs, + const bool& gradients) const { // Get pointer on CTA observation const GCTAObservation* ctaobs = dynamic_cast(&obs); if (ctaobs == NULL) { std::string msg = "Specified observation is not a CTA observation.\n" + obs.print(); - throw GException::invalid_argument(G_EVAL_GRADIENTS, msg); + throw GException::invalid_argument(G_EVAL, msg); } // Get pointer on CTA pointing @@ -384,11 +324,12 @@ double GCTAModelRadialAcceptance::eval_gradients(const GEvent& event, // Evaluate function and gradients double rad = (radial() != NULL) - ? radial()->eval_gradients(offset) : 1.0; + ? radial()->eval(offset, gradients) : 1.0; double spec = (spectral() != NULL) - ? spectral()->eval_gradients(event.energy(), event.time()) : 1.0; + ? spectral()->eval(event.energy(), event.time(), gradients) + : 1.0; double temp = (temporal() != NULL) - ? temporal()->eval_gradients(event.time()) : 1.0; + ? temporal()->eval(event.time(), gradients) : 1.0; // Compute value double value = rad * spec * temp; @@ -397,34 +338,39 @@ double GCTAModelRadialAcceptance::eval_gradients(const GEvent& event, double deadc = obs.deadc(event.time()); value *= deadc; - // Multiply factors to radial gradients - if (radial() != NULL) { - double fact = spec * temp * deadc; - if (fact != 1.0) { - for (int i = 0; i < radial()->size(); ++i) - (*radial())[i].factor_gradient( (*radial())[i].factor_gradient() * fact ); + // Optionally compute partial derivatives + if (gradients) { + + // Multiply factors to radial gradients + if (radial() != NULL) { + double fact = spec * temp * deadc; + if (fact != 1.0) { + for (int i = 0; i < radial()->size(); ++i) + (*radial())[i].factor_gradient((*radial())[i].factor_gradient() * fact ); + } } - } - // Multiply factors to spectral gradients - if (spectral() != NULL) { - double fact = rad * temp * deadc; - if (fact != 1.0) { - for (int i = 0; i < spectral()->size(); ++i) - (*spectral())[i].factor_gradient( (*spectral())[i].factor_gradient() * fact ); + // Multiply factors to spectral gradients + if (spectral() != NULL) { + double fact = rad * temp * deadc; + if (fact != 1.0) { + for (int i = 0; i < spectral()->size(); ++i) + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact ); + } } - } - // Multiply factors to temporal gradients - if (temporal() != NULL) { - double fact = rad * spec * deadc; - if (fact != 1.0) { - for (int i = 0; i < temporal()->size(); ++i) - (*temporal())[i].factor_gradient( (*temporal())[i].factor_gradient() * fact ); + // Multiply factors to temporal gradients + if (temporal() != NULL) { + double fact = rad * spec * deadc; + if (fact != 1.0) { + for (int i = 0; i < temporal()->size(); ++i) + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact ); + } } - } - // Return value + } // endif: computed partial derivatives + + // Return return value; } diff --git a/inst/cta/src/GCTAModelRadialGauss.cpp b/inst/cta/src/GCTAModelRadialGauss.cpp index 841c04f..3a7ff9b 100644 --- a/inst/cta/src/GCTAModelRadialGauss.cpp +++ b/inst/cta/src/GCTAModelRadialGauss.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialGauss.cpp - Radial Gaussian CTA model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -217,6 +217,8 @@ GCTAModelRadialGauss* GCTAModelRadialGauss::clone(void) const * @brief Evaluate function * * @param[in] offset Offset angle [degrees]. + * @param[in] gradients Compute gradients? + * @return Function value * * Evaluates the Gaussian model for a given offset. The Gaussian model is * defined as @@ -226,66 +228,36 @@ GCTAModelRadialGauss* GCTAModelRadialGauss::clone(void) const * \f$\theta\f$ is the offset angle (in degrees), and * \f$\sigma\f$ is the Gaussian width (in degrees\f$^2\f$). * - * Note that this method implements a function which is unity for - * \f$\theta=0\f$. - ***************************************************************************/ -double GCTAModelRadialGauss::eval(const double& offset) const -{ - // Compute value - double arg = offset * offset / sigma(); - double value = std::exp(-0.5 * arg * arg); - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GCTAModelRadialGauss::eval"; - std::cout << "(offset=" << offset << "): NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", sigma=" << sigma(); - std::cout << ")" << std::endl; - } - #endif - - // Return value - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] offset Offset angle [degrees]. - * - * Evaluates the Gaussian model for a given offset. The Gaussian model is - * defined as - * \f[f(\theta) = \exp \left(-\frac{1}{2} - * \left( \frac{\theta^2}{\sigma} \right)^2 \right)\f] - * where - * \f$\theta\f$ is the offset angle (in degrees), and - * \f$\sigma\f$ is the Gaussian width (in degrees\f$^2\f$). - * - * The partial derivative of the Gaussian width is given by + * If the @p gradients flag is true the method will also compute the partial + * derivatives of the parameters. The partial derivative of the Gaussian width + * is given by * \f[\frac{df}{d\sigma_v} = f(\theta) \frac{\theta^4}{\sigma^3} \sigma_s\f] * where * \f$\sigma_v\f$ is the value part, * \f$\sigma_s\f$ is the scaling part, and - * \f$\sigma = \sigma_v \sigma_s\f$. + * \f$\sigma = \sigma_v \sigma_s\f$. * * Note that this method implements a function which is unity for * \f$\theta=0\f$. ***************************************************************************/ -double GCTAModelRadialGauss::eval_gradients(const double& offset) const +double GCTAModelRadialGauss::eval(const double& offset, + const bool& gradients) const { // Compute value double arg = offset * offset / sigma(); double arg2 = arg * arg; - double value = exp(-0.5 * arg2); + double value = std::exp(-0.5 * arg2); + + // Optionally compute partial derivatives + if (gradients) { + + // Compute partial derivatives of the sigma parameter. + double g_sigma = value * arg2 / sigma() * m_sigma.scale(); - // Compute partial derivatives of the sigma parameter. - double g_sigma = value * arg2 / sigma() * m_sigma.scale(); + // Set factor gradient + m_sigma.factor_gradient(g_sigma); - // Set gradients (circumvent const correctness) - const_cast(this)->m_sigma.factor_gradient(g_sigma); + } // endif: computed partial derivatives // Compile option: Check for NaN/Inf #if defined(G_NAN_CHECK) diff --git a/inst/cta/src/GCTAModelRadialPolynom.cpp b/inst/cta/src/GCTAModelRadialPolynom.cpp index e114f87..8fabeec 100644 --- a/inst/cta/src/GCTAModelRadialPolynom.cpp +++ b/inst/cta/src/GCTAModelRadialPolynom.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialPolynom.cpp - Radial Polynom CTA model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -243,6 +243,8 @@ GCTAModelRadialPolynom* GCTAModelRadialPolynom::clone(void) const * @brief Evaluate function * * @param[in] offset Offset angle [degrees]. + * @param[in] gradients Compute gradients? + * @return Function value * * Evaluates the radial polynomial model for a given offset. The model is * defined as @@ -250,69 +252,18 @@ GCTAModelRadialPolynom* GCTAModelRadialPolynom::clone(void) const * where * \f$\theta\f$ is the offset angle (in degrees), and * \f$c_i\f$ are the polynomial coefficients. - ***************************************************************************/ -double GCTAModelRadialPolynom::eval(const double& offset) const -{ - // Initialise result - double value = 0.0; - - // Determine polynomial degree - int ncoeffs = m_coeffs.size(); - - // Compute value and gradients (only if there are coefficients) - if (ncoeffs > 0) { - - // Compute value - value = m_coeffs[ncoeffs-1].value(); - for (int i = ncoeffs-2; i >= 0; i--) { - value = value * offset + m_coeffs[i].value(); - } - - } // endif: there were coefficients - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GCTAModelRadialPolynom::eval"; - std::cout << "(offset=" << offset << "): NaN/Inf encountered"; - std::cout << " (value=" << value; - for (int i = 0; i < ncoeffs; ++i) { - std::cout << ", c" << i << "=" << m_coeffs[i].value(); - } - std::cout << ")" << std::endl; - } - #endif - - // Return value - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] offset Offset angle [degrees]. * - * Evaluates the Polynom model for a given offset. The Polynom model is - * defined as - * Evaluates the Polynom model for a given offset. The Polynomian model is - * defined as - * \f[f(\theta) = \sum_{i=0}^m c_i \theta^i\f] - * where - * \f$\theta\f$ is the offset angle (in degrees), and - * \f$c_i\f$ are the polynomial coefficients. - * - * The partial derivative of the Polynom are given by + * If the @p gradients flag is true the method will also compute the partial + * derivatives of the parameters. The partial derivative of the Polynom are + * given by * \f[\frac{df}{d{c_i}_v} = {c_i}_s \theta^i\f] * where * \f${c_i}_v\f$ is the value part, * \f${c_i}_s\f$ is the scaling part, and * \f$c_i = {c_i}_v {c_i}_s\f$. - * - * Note that this method implements a function which is unity for - * \f$\theta=0\f$. ***************************************************************************/ -double GCTAModelRadialPolynom::eval_gradients(const double& offset) const +double GCTAModelRadialPolynom::eval(const double& offset, + const bool& gradients) const { // Initialise result double value = 0.0; @@ -329,33 +280,38 @@ double GCTAModelRadialPolynom::eval_gradients(const double& offset) const value = value * offset + m_coeffs[i].value(); } - // Initialise theta^i for the first coefficient - double offset_power = 1.0; + // Optionally compute partial derivatives + if (gradients) { - // Compute gradients for all coefficients - for (int i = 0; i < ncoeffs; ++i) { + // Initialise theta^i for the first coefficient + double offset_power = 1.0; + + // Compute gradients for all coefficients + for (int i = 0; i < ncoeffs; ++i) { - // Compute gradient - double grad = offset_power * m_coeffs[i].scale(); + // Compute gradient + double grad = offset_power * m_coeffs[i].scale(); - // Store gradient - const_cast(this)->m_coeffs[i].factor_gradient(grad); + // Store gradient + m_coeffs[i].factor_gradient(grad); - // Increase offset power for next coefficient - offset_power *= offset; + // Increase offset power for next coefficient + offset_power *= offset; - } // endfor: looped over all coefficients + } // endfor: looped over all coefficients + + } // endif: computed partial derivatives } // endif: there were coefficients // Compile option: Check for NaN/Inf #if defined(G_NAN_CHECK) if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GCTAModelRadialPolynom::eval_gradients"; + std::cout << "*** ERROR: GCTAModelRadialPolynom::eval"; std::cout << "(offset=" << offset << "): NaN/Inf encountered"; std::cout << " (value=" << value; for (int i = 0; i < ncoeffs; ++i) { - std::cout << ", c_" << i << "=" << m_coeffs[i].value(); + std::cout << ", c" << i << "=" << m_coeffs[i].value(); } std::cout << ")" << std::endl; } diff --git a/inst/cta/src/GCTAModelRadialProfile.cpp b/inst/cta/src/GCTAModelRadialProfile.cpp index 648c4ee..003bf34 100644 --- a/inst/cta/src/GCTAModelRadialProfile.cpp +++ b/inst/cta/src/GCTAModelRadialProfile.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAModelRadialProfile.cpp - Radial Profile CTA model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -224,6 +224,8 @@ GCTAModelRadialProfile* GCTAModelRadialProfile::clone(void) const * @brief Evaluate function * * @param[in] offset Offset angle [degrees]. + * @param[in] gradients Compute gradients? + * @return Function value * * Evaluates the Profile model for a given offset. The Profile model is * defined as @@ -234,10 +236,11 @@ GCTAModelRadialProfile* GCTAModelRadialProfile::clone(void) const * \f$c_1\f$ is the width of the central plateau (core), and * \f$c_2\f$ is the size of the tail (tail). * - * Note that this method implements a function which is unity for - * \f$\theta=0\f$. + * Note that no analytical partial derivatives are implemented for this + * function. ***************************************************************************/ -double GCTAModelRadialProfile::eval(const double& offset) const +double GCTAModelRadialProfile::eval(const double& offset, + const bool& gradients) const { // Compute value double arg = 1.0 + std::pow(offset / width(), core()); @@ -262,22 +265,6 @@ double GCTAModelRadialProfile::eval(const double& offset) const /***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] offset Offset angle [degrees]. - * - * This method simply calls GCTAModelRadialProfile::eval() as no analytical - * gradients will be computed. See GCTAModelRadialProfile::eval() for details - * about the implemented method. - ***************************************************************************/ -double GCTAModelRadialProfile::eval_gradients(const double& offset) const -{ - // Return value - return (eval(offset)); -} - - -/***********************************************************************//** * @brief Returns MC instrument direction * * @param[in] dir Pointing direction. diff --git a/inst/cta/src/GCTAOnOffObservation.cpp b/inst/cta/src/GCTAOnOffObservation.cpp index 68621d7..c30f61c 100644 --- a/inst/cta/src/GCTAOnOffObservation.cpp +++ b/inst/cta/src/GCTAOnOffObservation.cpp @@ -1390,8 +1390,10 @@ double GCTAOnOffObservation::N_gamma(const GModels& models, // and the livetime (s) value += spectral->flux(emin, emax) * norm_flux; - // Determine the model gradients at the current energy - spectral->eval_gradients(emean); + // Determine the model gradients at the current energy. The + // eval() method needs a time in case that the spectral model + // has a time dependence. We simply use a dummy time here. + spectral->eval(emean, GTime(), true); // Loop over spectral model parameters for (int k = 0; k < spectral->size(); ++k, ++ipar) { @@ -1445,8 +1447,8 @@ double GCTAOnOffObservation::N_gamma(const GModels& models, * background models; that's not very satisfactory ... ***********************************************************************/ double GCTAOnOffObservation::N_bgd(const GModels& models, - const int& ibin, - GVector* grad) const + const int& ibin, + GVector* grad) const { // Get total number of model parameters int npars = models.npars(); @@ -1505,8 +1507,10 @@ double GCTAOnOffObservation::N_bgd(const GModels& models, // Determine the number of background events in model by // computing the model normalization at the mean bin energy bin // and multiplying the normalisation with the number of - // background events - value += spectral->eval_gradients(emean) * norm; + // background events. The eval() method needs a time in case + // that the spectral model has a time dependence. We simply + // use a dummy time here. + value += spectral->eval(emean, GTime(), true) * norm; // Compute the parameter gradients for all spectral model // parameters diff --git a/inst/cta/test/dev/test_gauss.py b/inst/cta/test/dev/test_gauss.py index bde91d3..a61d7f4 100755 --- a/inst/cta/test/dev/test_gauss.py +++ b/inst/cta/test/dev/test_gauss.py @@ -39,9 +39,9 @@ def gauss(dir, sigma): t = GTime() f1 = GModelSpatialRadialGauss(s, sigma) f2 = GModelSpatialRadialGauss(s, sigma + dh) - v1 = f1.eval_gradients(GPhoton(dir, e, t)) + v1 = f1.eval(GPhoton(dir, e, t), True) g1 = f1[2].gradient() - v2 = f2.eval_gradients(GPhoton(dir, e, t)) + v2 = f2.eval(GPhoton(dir, e, t), True) g2 = f2[2].gradient() g = (v2 - v1) / dh diff --git a/inst/cta/test/dev/test_radial_acceptance.py b/inst/cta/test/dev/test_radial_acceptance.py index b30df8f..8c381f7 100755 --- a/inst/cta/test/dev/test_radial_acceptance.py +++ b/inst/cta/test/dev/test_radial_acceptance.py @@ -39,7 +39,7 @@ def show_model(model): # Get model values values = [model.eval(theta) for theta in thetas] - values_grad = [model.eval_gradients(theta) for theta in thetas] + values_grad = [model.eval(theta, True) for theta in thetas] # Create figure plt.figure(1) diff --git a/pyext/GModel.i b/pyext/GModel.i index e7cce46..11881fb 100644 --- a/pyext/GModel.i +++ b/pyext/GModel.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModel.i - Abstract virtual model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2015 by Juergen Knoedlseder * + * copyright (C) 2009-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -51,9 +51,8 @@ public: virtual std::string type(void) const = 0; virtual bool is_constant(void) const = 0; virtual double eval(const GEvent& event, - const GObservation& obs) const = 0; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const = 0; + const GObservation& obs, + const bool& gradients = false) const = 0; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const = 0; virtual void read(const GXmlElement& xml) = 0; diff --git a/pyext/GModelData.i b/pyext/GModelData.i index e4d7169..c958f71 100644 --- a/pyext/GModelData.i +++ b/pyext/GModelData.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelData.i - Abstract virtual data model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -50,9 +50,8 @@ public: virtual std::string type(void) const = 0; virtual bool is_constant(void) const = 0; virtual double eval(const GEvent& event, - const GObservation& obs) const = 0; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const = 0; + const GObservation& obs, + const bool& gradients = false) const = 0; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const = 0; virtual GEvents* mc(const GObservation& obs, GRan& ran) const = 0; diff --git a/pyext/GModelSky.i b/pyext/GModelSky.i index 9ec6947..a7efcd5 100644 --- a/pyext/GModelSky.i +++ b/pyext/GModelSky.i @@ -61,9 +61,8 @@ public: virtual std::string type(void) const; virtual bool is_constant(void) const; virtual double eval(const GEvent& event, - const GObservation& obs) const; - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const; + const GObservation& obs, + const bool& gradients = false) const; virtual double npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const; diff --git a/pyext/GModelSpatial.i b/pyext/GModelSpatial.i index 2ebae0e..aafbdb3 100644 --- a/pyext/GModelSpatial.i +++ b/pyext/GModelSpatial.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatial.i - Spatial model abstract base class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2015 by Juergen Knoedlseder * + * copyright (C) 2009-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -47,8 +47,8 @@ public: virtual GModelSpatial* clone(void) const = 0; virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; - virtual double eval(const GPhoton& photon) const = 0; - virtual double eval_gradients(const GPhoton& photon) const = 0; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const = 0; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const = 0; virtual double mc_norm(const GSkyDir& dir, diff --git a/pyext/GModelSpatialDiffuse.i b/pyext/GModelSpatialDiffuse.i index b1e51ed..1884f72 100644 --- a/pyext/GModelSpatialDiffuse.i +++ b/pyext/GModelSpatialDiffuse.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialDiffuse.i - Abstract diffuse spatial model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2015 by Juergen Knoedlseder * + * copyright (C) 2013-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -46,8 +46,8 @@ public: virtual GModelSpatialDiffuse* clone(void) const = 0; virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; - virtual double eval(const GPhoton& photon) const = 0; - virtual double eval_gradients(const GPhoton& photon) const = 0; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const = 0; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const = 0; virtual double mc_norm(const GSkyDir& dir, diff --git a/pyext/GModelSpatialDiffuseConst.i b/pyext/GModelSpatialDiffuseConst.i index 2317e02..ec905d7 100644 --- a/pyext/GModelSpatialDiffuseConst.i +++ b/pyext/GModelSpatialDiffuseConst.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialDiffuseConst.i - Spatial isotropic model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -48,8 +48,8 @@ public: virtual GModelSpatialDiffuseConst* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/pyext/GModelSpatialDiffuseCube.i b/pyext/GModelSpatialDiffuseCube.i index 489da83..3e5f07b 100644 --- a/pyext/GModelSpatialDiffuseCube.i +++ b/pyext/GModelSpatialDiffuseCube.i @@ -53,11 +53,11 @@ public: virtual GModelSpatialDiffuseCube* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, - const GTime& time, - GRan& ran) const; + const GTime& time, + GRan& ran) const; virtual double mc_norm(const GSkyDir& dir, const double& radius) const; virtual bool contains(const GSkyDir& dir, diff --git a/pyext/GModelSpatialDiffuseMap.i b/pyext/GModelSpatialDiffuseMap.i index f44f834..c41a23f 100644 --- a/pyext/GModelSpatialDiffuseMap.i +++ b/pyext/GModelSpatialDiffuseMap.i @@ -53,8 +53,8 @@ public: virtual GModelSpatialDiffuseMap* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/pyext/GModelSpatialElliptical.i b/pyext/GModelSpatialElliptical.i index 3a282fd..fc297f4 100644 --- a/pyext/GModelSpatialElliptical.i +++ b/pyext/GModelSpatialElliptical.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialElliptical.i - Abstract elliptical spatial model class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2015 by Juergen Knoedlseder * + * copyright (C) 2013-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -50,21 +50,18 @@ public: virtual double eval(const double& theta, const double& posangle, const GEnergy& energy, - const GTime& time) const = 0; - virtual double eval_gradients(const double& theta, - const double& posangle, - const GEnergy& energy, - const GTime& time) const = 0; + const GTime& time, + const bool& gradients = false) const = 0; virtual bool contains(const GSkyDir& dir, const double& margin = 0.0) const = 0; virtual GSkyDir mc(const GEnergy& energy, - const GTime& time, - GRan& ran) const = 0; + const GTime& time, + GRan& ran) const = 0; virtual double theta_max(void) const = 0; // Implemented virtual methods - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual double mc_norm(const GSkyDir& dir, const double& radius) const; virtual void read(const GXmlElement& xml); virtual void write(GXmlElement& xml) const; diff --git a/pyext/GModelSpatialEllipticalDisk.i b/pyext/GModelSpatialEllipticalDisk.i index 1900664..68fa402 100644 --- a/pyext/GModelSpatialEllipticalDisk.i +++ b/pyext/GModelSpatialEllipticalDisk.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialEllipticalDisk.i - Elliptical disk source model class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2015 by Michael Mayer * + * copyright (C) 2013-2016 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -55,11 +55,8 @@ public: virtual double eval(const double& theta, const double& posangle, const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const double& posangle, - const GEnergy& energy, - const GTime& time) const; + const GTime& time, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; @@ -78,10 +75,7 @@ public: GModelSpatialEllipticalDisk copy() { return (*self); } - double eval(const GPhoton& photon) const { - return self->GModelSpatialElliptical::eval(photon); - } - double eval_gradients(const GPhoton& photon) const { - return self->GModelSpatialElliptical::eval_gradients(photon); + double eval(const GPhoton& photon, const bool& gradients) const { + return self->GModelSpatialElliptical::eval(photon, gradients); } }; diff --git a/pyext/GModelSpatialEllipticalGauss.i b/pyext/GModelSpatialEllipticalGauss.i index ac9508f..222ee0c 100644 --- a/pyext/GModelSpatialEllipticalGauss.i +++ b/pyext/GModelSpatialEllipticalGauss.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialEllipticalGauss.i - Elliptical gauss source model class * * ----------------------------------------------------------------------- * - * copyright (C) 2015 by Michael Mayer * + * copyright (C) 2015-2016 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -55,11 +55,8 @@ public: virtual double eval(const double& theta, const double& posangle, const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const double& posangle, - const GEnergy& energy, - const GTime& time) const; + const GTime& time, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; @@ -78,10 +75,7 @@ public: GModelSpatialEllipticalGauss copy() { return (*self); } - double eval(const GPhoton& photon) const { - return self->GModelSpatialElliptical::eval(photon); - } - double eval_gradients(const GPhoton& photon) const { - return self->GModelSpatialElliptical::eval_gradients(photon); + double eval(const GPhoton& photon, const bool& gradients) const { + return self->GModelSpatialElliptical::eval(photon, gradients); } }; diff --git a/pyext/GModelSpatialPointSource.i b/pyext/GModelSpatialPointSource.i index c91e091..97a0654 100644 --- a/pyext/GModelSpatialPointSource.i +++ b/pyext/GModelSpatialPointSource.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialPointSource.i - Spatial point source model class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2015 by Juergen Knoedlseder * + * copyright (C) 2009-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -50,8 +50,8 @@ public: virtual GModelSpatialPointSource* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; diff --git a/pyext/GModelSpatialRadial.i b/pyext/GModelSpatialRadial.i index 6962710..887b604 100644 --- a/pyext/GModelSpatialRadial.i +++ b/pyext/GModelSpatialRadial.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialRadial.i - Abstract radial spatial model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -50,10 +50,8 @@ public: virtual std::string type(void) const = 0; virtual double eval(const double& theta, const GEnergy& energy, - const GTime& time) const = 0; - virtual double eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const = 0; + const GTime& time, + const bool& gradients = false) const = 0; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const = 0; @@ -62,8 +60,8 @@ public: virtual double theta_max(void) const = 0; // Implemented virtual base class methods - virtual double eval(const GPhoton& photon) const; - virtual double eval_gradients(const GPhoton& photon) const; + virtual double eval(const GPhoton& photon, + const bool& gradients = false) const; virtual double mc_norm(const GSkyDir& dir, const double& radius) const; virtual void read(const GXmlElement& xml); virtual void write(GXmlElement& xml) const; diff --git a/pyext/GModelSpatialRadialDisk.i b/pyext/GModelSpatialRadialDisk.i index 551fb25..9f803cd 100644 --- a/pyext/GModelSpatialRadialDisk.i +++ b/pyext/GModelSpatialRadialDisk.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialRadialDisk.i - Radial disk source model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Christoph Deil * + * copyright (C) 2011-2016 by Christoph Deil * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -54,10 +54,8 @@ public: virtual std::string type(void) const; virtual double eval(const double& theta, const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const; + const GTime& time, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; @@ -76,9 +74,9 @@ public: /***********************************************************************//** * @brief GModelSpatialRadialDisk class extension * - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the - * extension to force swig to build also the interface for these methods that - * are implemented in the base class only. It's not clear to me why these + * The eval(GPhoton&, bool&) method need to be defined in the extension to + * force swig to build also the interface for these methods that are + * implemented in the base class only. It's not clear to me why these * methods are not inherited automatically. Maybe this could also be handled * by a %typemap(typecheck) construct. ***************************************************************************/ @@ -86,10 +84,7 @@ public: GModelSpatialRadialDisk copy() { return (*self); } - double eval(const GPhoton& photon) const { - return self->GModelSpatialRadial::eval(photon); - } - double eval_gradients(const GPhoton& photon) const { - return self->GModelSpatialRadial::eval_gradients(photon); + double eval(const GPhoton& photon, const bool& gradients) const { + return self->GModelSpatialRadial::eval(photon, gradients); } }; diff --git a/pyext/GModelSpatialRadialGauss.i b/pyext/GModelSpatialRadialGauss.i index 06d70e5..71cb843 100644 --- a/pyext/GModelSpatialRadialGauss.i +++ b/pyext/GModelSpatialRadialGauss.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialRadialGauss.i - Radial Gaussian source model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -51,10 +51,8 @@ public: virtual std::string type(void) const; virtual double eval(const double& theta, const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const; + const GTime& time, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; @@ -73,9 +71,9 @@ public: /***********************************************************************//** * @brief GModelSpatialRadialGauss class extension * - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the - * extension to force swig to build also the interface for these methods that - * are implemented in the base class only. It's not clear to me why these + * The eval(GPhoton&, bool&) method need to be defined in the extension to + * force swig to build also the interface for these methods that are + * implemented in the base class only. It's not clear to me why these * methods are not inherited automatically. Maybe this could also be handled * by a %typemap(typecheck) construct. ***************************************************************************/ @@ -83,10 +81,7 @@ public: GModelSpatialRadialGauss copy() { return (*self); } - double eval(const GPhoton& photon) const { - return self->GModelSpatialRadial::eval(photon); - } - double eval_gradients(const GPhoton& photon) const { - return self->GModelSpatialRadial::eval_gradients(photon); + double eval(const GPhoton& photon, const bool& gradients) const { + return self->GModelSpatialRadial::eval(photon, gradients); } }; diff --git a/pyext/GModelSpatialRadialShell.i b/pyext/GModelSpatialRadialShell.i index 51db8fa..5cab68c 100644 --- a/pyext/GModelSpatialRadialShell.i +++ b/pyext/GModelSpatialRadialShell.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialRadialShell.i - Radial spatial shell source model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2015 by Christoph Deil * + * copyright (C) 2011-2016 by Christoph Deil * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -52,10 +52,8 @@ public: virtual std::string type(void) const; virtual double eval(const double& theta, const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const; + const GTime& time, + const bool& gradients = false) const; virtual GSkyDir mc(const GEnergy& energy, const GTime& time, GRan& ran) const; @@ -76,9 +74,9 @@ public: /***********************************************************************//** * @brief GModelSpatialRadialShell class extension * - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the - * extension to force swig to build also the interface for these methods that - * are implemented in the base class only. It's not clear to me why these + * The eval(GPhoton&, bool&) method need to be defined in the extension to + * force swig to build also the interface for these methods that are + * implemented in the base class only. It's not clear to me why these * methods are not inherited automatically. Maybe this could also be handled * by a %typemap(typecheck) construct. ***************************************************************************/ @@ -86,10 +84,7 @@ public: GModelSpatialRadialShell copy() { return (*self); } - double eval(const GPhoton& photon) const { - return self->GModelSpatialRadial::eval(photon); - } - double eval_gradients(const GPhoton& photon) const { - return self->GModelSpatialRadial::eval_gradients(photon); + double eval(const GPhoton& photon, const bool& gradients) const { + return self->GModelSpatialRadial::eval(photon, gradients); } }; diff --git a/pyext/GModelSpectral.i b/pyext/GModelSpectral.i index b7d262f..6aa133e 100644 --- a/pyext/GModelSpectral.i +++ b/pyext/GModelSpectral.i @@ -48,9 +48,8 @@ public: virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const = 0; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()) = 0; + const GTime& srcTime = GTime(), + const bool& gradients = false) const = 0; virtual double flux(const GEnergy& emin, const GEnergy& emax) const = 0; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralBrokenPlaw.i b/pyext/GModelSpectralBrokenPlaw.i index fe2666d..2ffa8fa 100644 --- a/pyext/GModelSpectralBrokenPlaw.i +++ b/pyext/GModelSpectralBrokenPlaw.i @@ -57,9 +57,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralConst.i b/pyext/GModelSpectralConst.i index 2bb9a86..0c8ae8c 100644 --- a/pyext/GModelSpectralConst.i +++ b/pyext/GModelSpectralConst.i @@ -50,9 +50,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralExpInvPlaw.i b/pyext/GModelSpectralExpInvPlaw.i index 8a7b983..9fb37fc 100644 --- a/pyext/GModelSpectralExpInvPlaw.i +++ b/pyext/GModelSpectralExpInvPlaw.i @@ -56,18 +56,14 @@ public: GModelSpectralExpInvPlaw(const GModelSpectralExpInvPlaw& model); virtual ~GModelSpectralExpInvPlaw(void); - // Operators - virtual GModelSpectralExpInvPlaw& operator=(const GModelSpectralExpInvPlaw& model); - // Implemented pure virtual methods virtual void clear(void); virtual GModelSpectralExpInvPlaw* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralExpPlaw.i b/pyext/GModelSpectralExpPlaw.i index 8005217..903b497 100644 --- a/pyext/GModelSpectralExpPlaw.i +++ b/pyext/GModelSpectralExpPlaw.i @@ -57,9 +57,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralFunc.i b/pyext/GModelSpectralFunc.i index 66e872f..a059197 100644 --- a/pyext/GModelSpectralFunc.i +++ b/pyext/GModelSpectralFunc.i @@ -50,9 +50,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralGauss.i b/pyext/GModelSpectralGauss.i index 4162c95..8752e29 100644 --- a/pyext/GModelSpectralGauss.i +++ b/pyext/GModelSpectralGauss.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpectralGauss.i - Spectral gaussian model class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2015 by Christoph Deil & Ellis Owen * + * copyright (C) 2014-2016 by Christoph Deil & Ellis Owen * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -51,9 +51,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralLogParabola.i b/pyext/GModelSpectralLogParabola.i index 9632155..a511aeb 100644 --- a/pyext/GModelSpectralLogParabola.i +++ b/pyext/GModelSpectralLogParabola.i @@ -57,9 +57,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralNodes.i b/pyext/GModelSpectralNodes.i index 91f6077..a00aacd 100644 --- a/pyext/GModelSpectralNodes.i +++ b/pyext/GModelSpectralNodes.i @@ -48,9 +48,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralPlaw.i b/pyext/GModelSpectralPlaw.i index 5ee7fff..4e1c3c9 100644 --- a/pyext/GModelSpectralPlaw.i +++ b/pyext/GModelSpectralPlaw.i @@ -55,9 +55,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralPlawEnergyFlux.i b/pyext/GModelSpectralPlawEnergyFlux.i index fb7b492..59b2c24 100644 --- a/pyext/GModelSpectralPlawEnergyFlux.i +++ b/pyext/GModelSpectralPlawEnergyFlux.i @@ -57,9 +57,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralPlawPhotonFlux.i b/pyext/GModelSpectralPlawPhotonFlux.i index 103fd18..ff2b2fb 100644 --- a/pyext/GModelSpectralPlawPhotonFlux.i +++ b/pyext/GModelSpectralPlawPhotonFlux.i @@ -57,9 +57,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelSpectralSuperExpPlaw.i b/pyext/GModelSpectralSuperExpPlaw.i index 8018c9b..2458292 100644 --- a/pyext/GModelSpectralSuperExpPlaw.i +++ b/pyext/GModelSpectralSuperExpPlaw.i @@ -59,9 +59,8 @@ public: virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); + const GTime& srcTime = GTime(), + const bool& gradients = false) const; virtual double flux(const GEnergy& emin, const GEnergy& emax) const; virtual double eflux(const GEnergy& emin, diff --git a/pyext/GModelTemporal.i b/pyext/GModelTemporal.i index 60d1ab0..76c6780 100644 --- a/pyext/GModelTemporal.i +++ b/pyext/GModelTemporal.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelTemporal.i - Abstract temporal model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2015 by Juergen Knoedlseder * + * copyright (C) 2009-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -47,10 +47,12 @@ public: virtual GModelTemporal* clone(void) const = 0; virtual std::string classname(void) const = 0; virtual std::string type(void) const = 0; - virtual double eval(const GTime& srcTime) const = 0; - virtual double eval_gradients(const GTime& srcTime) = 0; - virtual GTimes mc(const double& rate, const GTime& tmin, - const GTime& tmax, GRan& ran) const = 0; + virtual double eval(const GTime& srcTime, + const bool& gradients = false) const = 0; + virtual GTimes mc(const double& rate, + const GTime& tmin, + const GTime& tmax, + GRan& ran) const = 0; virtual void read(const GXmlElement& xml) = 0; virtual void write(GXmlElement& xml) const = 0; diff --git a/pyext/GModelTemporalConst.i b/pyext/GModelTemporalConst.i index 48b0b62..6a84e8c 100644 --- a/pyext/GModelTemporalConst.i +++ b/pyext/GModelTemporalConst.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelTemporalConst.i - Temporal constant model class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2015 by Juergen Knoedlseder * + * copyright (C) 2009-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -48,10 +48,12 @@ public: virtual GModelTemporalConst* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; - virtual double eval(const GTime& srcTime) const; - virtual double eval_gradients(const GTime& srcTime); - virtual GTimes mc(const double& rate, const GTime& tmin, - const GTime& tmax, GRan& ran) const; + virtual double eval(const GTime& srcTime, + const bool& gradients = false) const; + virtual GTimes mc(const double& rate, + const GTime& tmin, + const GTime& tmax, + GRan& ran) const; virtual void read(const GXmlElement& xml); virtual void write(GXmlElement& xml) const; diff --git a/pyext/GModels.i b/pyext/GModels.i index f5652d1..05c0597 100644 --- a/pyext/GModels.i +++ b/pyext/GModels.i @@ -68,8 +68,9 @@ public: void write(GXml& xml) const; int npars(void) const; GOptimizerPars pars(void) const; - double eval(const GEvent& event, const GObservation& obs) const; - double eval_gradients(const GEvent& event, const GObservation& obs) const; + double eval(const GEvent& event, + const GObservation& obs, + const bool& gradients = false) const; }; diff --git a/src/model/GModelSky.cpp b/src/model/GModelSky.cpp index 3276e60..5dcc2b3 100644 --- a/src/model/GModelSky.cpp +++ b/src/model/GModelSky.cpp @@ -516,10 +516,10 @@ double GModelSky::value(const GPhoton& photon) GVector GModelSky::gradients(const GPhoton& photon) { // Evaluate source model gradients - if (m_spatial != NULL) m_spatial->eval_gradients(photon); - if (m_spectral != NULL) m_spectral->eval_gradients(photon.energy(), - photon.time()); - if (m_temporal != NULL) m_temporal->eval_gradients(photon.time()); + if (m_spatial != NULL) m_spatial->eval(photon, true); + if (m_spectral != NULL) m_spectral->eval(photon.energy(), photon.time(), + true); + if (m_temporal != NULL) m_temporal->eval(photon.time(), true); // Set vector of gradients GVector gradients; @@ -540,40 +540,21 @@ GVector GModelSky::gradients(const GPhoton& photon) * * @param[in] event Observed event. * @param[in] obs Observation. - * @return Value of sky model + * @param[in] gradients Compute gradients? + * @return Value of sky model. * * Evalutes the value of the sky model for an @p event of a specific * observation @p obs. - ***************************************************************************/ -double GModelSky::eval(const GEvent& event, const GObservation& obs) const -{ - // Evaluate function - double value = obs.response()->convolve(*this, event, obs, false); - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate sky model and parameter gradients for a given event of an - * observation - * - * @param[in] event Observed event. - * @param[in] obs Observation. - * @return Value of sky model - * - * Evalutes the value of the sky model and of the parameter for an @p event - * of a specific observation @p obs. * - * While the value of the sky model is returned by the method, the parameter - * gradients are set as GModelPar members. + * If the @p gradients flag is true the method will also compute the + * parameter gradients for all model parameters. ***************************************************************************/ -double GModelSky::eval_gradients(const GEvent& event, - const GObservation& obs) const +double GModelSky::eval(const GEvent& event, + const GObservation& obs, + const bool& gradients) const { // Evaluate function - double value = obs.response()->convolve(*this, event, obs, true); + double value = obs.response()->convolve(*this, event, obs, gradients); // Return return value; diff --git a/src/model/GModelSpatialDiffuseConst.cpp b/src/model/GModelSpatialDiffuseConst.cpp index b95efaa..d7246ae 100644 --- a/src/model/GModelSpatialDiffuseConst.cpp +++ b/src/model/GModelSpatialDiffuseConst.cpp @@ -272,6 +272,7 @@ GModelSpatialDiffuseConst* GModelSpatialDiffuseConst::clone(void) const * @brief Evaluate isotropic spatial model value * * @param[in] photon Incident photon. + * @param[in] gradients Compute gradients? * @return Model value (\f$sr^{-1}\f$). * * Evaluates the isotropic spatial model for a given @p photon, characterised @@ -285,39 +286,9 @@ GModelSpatialDiffuseConst* GModelSpatialDiffuseConst::clone(void) const * * where \f$v\f$ is the value() parameter, which is divided by the solid * angle \f$4\pi\f$ of the celestial sphere. - ***************************************************************************/ -double GModelSpatialDiffuseConst::eval(const GPhoton& photon) const -{ - // Set normalization constant - const double norm = 1.0 / gammalib::fourpi; - - // Compute model value - double value = norm * m_value.value(); - - // Return model value - return (value); -} - - -/***********************************************************************//** - * @brief Evaluate isotropic spatial model value and parameter gradient * - * @param[in] photon Incident photon. - * @return Model value (\f$sr^{-1}\f$). - * - * Evaluates the isotropic spatial model for a given @p photon, characterised - * by a sky direction, energy and arrival time \f$(\vec{p},E,t)\f$, and - * computes the gradient of the value() parameter. By definition, the - * isotropic spatial model and gradient are independent of sky direction, - * energy and time. The model value is given by - * - * \f[ - * M_\mathrm{S}(\vec{p}|E,t) = \frac{v}{4\pi} - * \f] - * - * where \f$v\f$ is the value() parameter, which is divided by the solid - * angle \f$4\pi\f$ of the celestial sphere. The value() gradient is given - * by + * If the @p gradients flag is true the method will also evaluate the partial + * derivatives of the model. The value() gradient is given by * * \f[ * \frac{\partial M_\mathrm{S}(\vec{p}|E,t)}{\partial v_v} = @@ -327,7 +298,8 @@ double GModelSpatialDiffuseConst::eval(const GPhoton& photon) const * where \f$v=v_v v_s\f$ is the factorisation of the value() parameter into * a value \f$v_v\f$ and scale \f$v_s\f$ term. ***************************************************************************/ -double GModelSpatialDiffuseConst::eval_gradients(const GPhoton& photon) const +double GModelSpatialDiffuseConst::eval(const GPhoton& photon, + const bool& gradients) const { // Set normalization constant const double norm = 1.0 / gammalib::fourpi; @@ -335,14 +307,19 @@ double GModelSpatialDiffuseConst::eval_gradients(const GPhoton& photon) const // Compute model value double value = norm * m_value.value(); - // Compute partial derivatives of the parameter values - double g_norm = (m_value.is_free()) ? norm * m_value.scale() : 0.0; + // Optionally compute partial derivatives + if (gradients) { + + // Compute partial derivatives of the parameter values + double g_norm = (m_value.is_free()) ? norm * m_value.scale() : 0.0; - // Set gradient (circumvent const correctness) - const_cast(this)->m_value.factor_gradient(g_norm); + // Set gradient + m_value.factor_gradient(g_norm); - // Return value - return value; + } // endif: computed partial derivatives + + // Return model value + return (value); } diff --git a/src/model/GModelSpatialDiffuseCube.cpp b/src/model/GModelSpatialDiffuseCube.cpp index c63691c..850eab5 100644 --- a/src/model/GModelSpatialDiffuseCube.cpp +++ b/src/model/GModelSpatialDiffuseCube.cpp @@ -310,11 +310,16 @@ GModelSpatialDiffuseCube* GModelSpatialDiffuseCube::clone(void) const * @brief Evaluate function * * @param[in] photon Incident photon. - * @return Model value. + * @param[in] gradients Compute gradients? + * @return Sky map intensity (\f$\mbox{ph cm}^{-2}\mbox{sr}^{-1}\mbox{s}^{-1}\f$) * * Computes the spatial diffuse model as function of photon parameters. + * + * If the @p gradients flag is true the method will also evaluate the partial + * derivatives of the model. ***************************************************************************/ -double GModelSpatialDiffuseCube::eval(const GPhoton& photon) const +double GModelSpatialDiffuseCube::eval(const GPhoton& photon, + const bool& gradients) const { // Get log-log interpolated cube intensity double intensity = cube_intensity(photon); @@ -322,45 +327,26 @@ double GModelSpatialDiffuseCube::eval(const GPhoton& photon) const // Set the intensity times the scaling factor as model value double value = intensity * m_value.value(); - // Make sure that value is not negative - if (value < 0.0) { - value = 0.0; - } + // Optionally compute partial derivatives + if (gradients) { - // Return value - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] photon Incident photon. - * @return Model value. - * - * Computes the spatial diffuse model as function of photon parameters and - * sets the value gradient. - ***************************************************************************/ -double GModelSpatialDiffuseCube::eval_gradients(const GPhoton& photon) const -{ - // Get log-log interpolated cube intensity - double intensity = cube_intensity(photon); + // Compute partial derivatives of the parameter value. In case that + // the value is negative set the gradient to zero. + double g_value = (m_value.is_free()) ? intensity * m_value.scale() : 0.0; + if (value < 0.0) { + g_value = 0.0; + } - // Compute the model value - double value = intensity * m_value.value(); + // Set gradient + m_value.factor_gradient(g_value); - // Compute partial derivatives of the parameter value - double g_value = (m_value.is_free()) ? intensity * m_value.scale() : 0.0; + } // endif: computed partial derivatives // Make sure that value is not negative if (value < 0.0) { - value = 0.0; - g_value = 0.0; + value = 0.0; } - // Set gradient to 0 (circumvent const correctness) - const_cast(this)->m_value.factor_gradient(g_value); - // Return value return value; } diff --git a/src/model/GModelSpatialDiffuseMap.cpp b/src/model/GModelSpatialDiffuseMap.cpp index e612d79..c6943d3 100644 --- a/src/model/GModelSpatialDiffuseMap.cpp +++ b/src/model/GModelSpatialDiffuseMap.cpp @@ -309,43 +309,32 @@ GModelSpatialDiffuseMap* GModelSpatialDiffuseMap::clone(void) const * @brief Return intensity of sky map * * @param[in] photon Incident photon. + * @param[in] gradients Compute gradients? * @return Sky map intensity (\f$\mbox{ph cm}^{-2}\mbox{sr}^{-1}\mbox{s}^{-1}\f$) * * Returns the intensity of the sky map at the specified sky direction * multiplied by the normalization factor. If the sky direction falls outside * the sky map, an intensity of 0 is returned. + * + * If the @p gradients flag is true the method will also evaluate the partial + * derivatives of the model. ***************************************************************************/ -double GModelSpatialDiffuseMap::eval(const GPhoton& photon) const +double GModelSpatialDiffuseMap::eval(const GPhoton& photon, + const bool& gradients) const { // Get skymap intensity double intensity = m_map(photon.dir()); - // Return intensity times normalization factor - return (intensity * m_value.value()); -} - + // Optionally compute partial derivatives + if (gradients) { -/***********************************************************************//** - * @brief Return intensity of sky map and gradient - * - * @param[in] photon Incident photon. - * @return Sky map intensity (\f$\mbox{ph cm}^{-2}\mbox{sr}^{-1}\mbox{s}^{-1}\f$) - * - * Returns the intensity of the sky map at the specified sky direction - * multiplied by the normalization factor. The method also sets the gradient - * with respect to the normalization factor. If the sky direction falls - * outside the sky map, an intensity of 0 is returned. - ***************************************************************************/ -double GModelSpatialDiffuseMap::eval_gradients(const GPhoton& photon) const -{ - // Get skymap intensity - double intensity = m_map(photon.dir()); + // Compute partial derivatives of the parameter values + double g_value = (m_value.is_free()) ? intensity * m_value.scale() : 0.0; - // Compute partial derivatives of the parameter values - double g_value = (m_value.is_free()) ? intensity * m_value.scale() : 0.0; + // Set gradient + m_value.factor_gradient(g_value); - // Set gradient to 0 (circumvent const correctness) - const_cast(this)->m_value.factor_gradient(g_value); + } // endif: computed partial derivatives // Return intensity times normalization factor return (intensity * m_value.value()); diff --git a/src/model/GModelSpatialElliptical.cpp b/src/model/GModelSpatialElliptical.cpp index c86f9d3..ff02fba 100644 --- a/src/model/GModelSpatialElliptical.cpp +++ b/src/model/GModelSpatialElliptical.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialElliptical.cpp - Abstract elliptical spatial model class * * ----------------------------------------------------------------------- * - * copyright (C) 2013-2014 by Juergen Knoedlseder * + * copyright (C) 2013-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -161,43 +161,25 @@ GModelSpatialElliptical& GModelSpatialElliptical::operator=(const GModelSpatialE * @brief Return model value * * @param[in] photon Incident Photon. + * @param[in] gradients Compute gradients? + * @return Value of spatial elliptical model. * * Evaluates the elliptical spatial model value for a specific incident * @p photon. - ***************************************************************************/ -double GModelSpatialElliptical::eval(const GPhoton& photon) const -{ - // Compute distance from source and position angle (in radians) - const GSkyDir& srcDir = photon.dir(); - double theta = dir().dist(srcDir); - double posang = dir().posang(srcDir); - - // Evaluate model - double value = eval(theta, posang, photon.energy(), photon.time()); - - // Return result - return value; -} - - -/***********************************************************************//** - * @brief Return model value and set analytical gradients - * - * @param[in] photon Incident Photon. * - * Evaluates the elliptical spatial model value and analytical model - * parameter gradients for a specific incident @p photon. + * If the @p gradients flag is true the method will also compute the + * parameter gradients for all model parameters. ***************************************************************************/ -double GModelSpatialElliptical::eval_gradients(const GPhoton& photon) const +double GModelSpatialElliptical::eval(const GPhoton& photon, + const bool& gradients) const { // Compute distance from source and position angle (in radians) const GSkyDir& srcDir = photon.dir(); - double theta = dir().dist(srcDir); - double posang = dir().posang(srcDir); + double theta = dir().dist(srcDir); + double posang = dir().posang(srcDir); - // Evaluate model and set gradients - double value = eval_gradients(theta, posang, photon.energy(), - photon.time()); + // Evaluate model + double value = eval(theta, posang, photon.energy(), photon.time(), gradients); // Return result return value; diff --git a/src/model/GModelSpatialEllipticalDisk.cpp b/src/model/GModelSpatialEllipticalDisk.cpp index b009706..79fa85f 100644 --- a/src/model/GModelSpatialEllipticalDisk.cpp +++ b/src/model/GModelSpatialEllipticalDisk.cpp @@ -240,6 +240,7 @@ GModelSpatialEllipticalDisk* GModelSpatialEllipticalDisk::clone(void) const * @param[in] posangle Position angle (counterclockwise from North) (radians). * @param[in] energy Photon energy. * @param[in] time Photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value. * * Evaluates the spatial component for an elliptical disk source model. The @@ -277,6 +278,10 @@ GModelSpatialEllipticalDisk* GModelSpatialEllipticalDisk::clone(void) const * The normalisation constant \f${\tt m\_norm}\f$ which is the inverse of the * solid angle subtended by an ellipse is given by * + * The method will not compute analytical parameter gradients, even if the + * @p gradients argument is set to true. Radial disk parameter gradients + * need to be computed numerically. + * * @todo Quote formula for ellipse solid angle * * (see the update() method). @@ -284,7 +289,8 @@ GModelSpatialEllipticalDisk* GModelSpatialEllipticalDisk::clone(void) const double GModelSpatialEllipticalDisk::eval(const double& theta, const double& posangle, const GEnergy& energy, - const GTime& time) const + const GTime& time, + const bool& gradients) const { // Initialise value double value = 0.0; @@ -329,31 +335,6 @@ double GModelSpatialEllipticalDisk::eval(const double& theta, /***********************************************************************//** - * @brief Evaluate function and gradients (in units of sr^-1) - * - * @param[in] theta Angular distance from disk centre (radians). - * @param[in] posangle Position angle (counterclockwise from North) (radians). - * @param[in] energy Photon energy. - * @param[in] time Photon arrival time. - * @return Model value. - * - * Evaluates the function value. No gradient computation is implemented as - * Elliptical models will be convolved with the instrument response and thus - * require the numerical computation of the derivatives. - * - * See the eval() method for more information. - ***************************************************************************/ -double GModelSpatialEllipticalDisk::eval_gradients(const double& theta, - const double& posangle, - const GEnergy& energy, - const GTime& time) const -{ - // Return value - return (eval(theta, posangle, energy, time)); -} - - -/***********************************************************************//** * @brief Returns MC sky direction * * @param[in] energy Photon energy. diff --git a/src/model/GModelSpatialEllipticalGauss.cpp b/src/model/GModelSpatialEllipticalGauss.cpp index 99236bb..02c5cbe 100644 --- a/src/model/GModelSpatialEllipticalGauss.cpp +++ b/src/model/GModelSpatialEllipticalGauss.cpp @@ -273,6 +273,7 @@ GModelSpatialEllipticalGauss* GModelSpatialEllipticalGauss::clone(void) const * @param[in] posangle Position angle (counterclockwise from North) (radians). * @param[in] energy Photon energy. * @param[in] time Photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value. * * Evaluates the spatial component for an elliptical Gaussian source model. @@ -309,6 +310,10 @@ GModelSpatialEllipticalGauss* GModelSpatialEllipticalGauss::clone(void) const * {\tt m\_norm} = \frac{1}{2 \pi \times a \times b} * \f] * + * The method will not compute analytical parameter gradients, even if the + * @p gradients argument is set to true. Radial disk parameter gradients + * need to be computed numerically. + * * @warning * The above formula for an elliptical Gaussian are accurate for small * angles, with semimajor and semiminor axes below a few degrees. This @@ -323,7 +328,8 @@ GModelSpatialEllipticalGauss* GModelSpatialEllipticalGauss::clone(void) const double GModelSpatialEllipticalGauss::eval(const double& theta, const double& posangle, const GEnergy& energy, - const GTime& time) const + const GTime& time, + const bool& gradients) const { // Initialise value double value = 0.0; @@ -386,31 +392,6 @@ double GModelSpatialEllipticalGauss::eval(const double& theta, /***********************************************************************//** - * @brief Evaluate function and gradients (in units of sr^-1) - * - * @param[in] theta Angular distance from gaussian centre (radians). - * @param[in] posangle Position angle (counterclockwise from North) (radians). - * @param[in] energy Photon energy. - * @param[in] time Photon arrival time. - * @return Model value. - * - * Evaluates the function value. No gradient computation is implemented as - * elliptical models will be convolved with the instrument response and thus - * require the numerical computation of the derivatives. - * - * See the eval() method for more information. - ***************************************************************************/ -double GModelSpatialEllipticalGauss::eval_gradients(const double& theta, - const double& posangle, - const GEnergy& energy, - const GTime& time) const -{ - // Return value - return (eval(theta, posangle, energy, time)); -} - - -/***********************************************************************//** * @brief Returns MC sky direction * * @param[in] energy Photon energy. diff --git a/src/model/GModelSpatialPointSource.cpp b/src/model/GModelSpatialPointSource.cpp index a59dbc1..0c33943 100644 --- a/src/model/GModelSpatialPointSource.cpp +++ b/src/model/GModelSpatialPointSource.cpp @@ -278,14 +278,20 @@ GModelSpatialPointSource* GModelSpatialPointSource::clone(void) const * @brief Evaluate function * * @param[in] photon Incident photon. + * @param[in] gradients Compute gradients? * @return Model value. * * Evaluates the spatial part for a point source model. It implements a delta * function with respect to the coordinates of the source. For numerical * reasons, a certain tolerance is accepted (typically 0.1 arcsec, i.e. well * below the angular resolution of gamma-ray telescopes). + * + * The method will not compute analytical parameter gradients, even if the + * @p gradients argument is set to true. Point source parameter gradients + * need to be computed numerically. ***************************************************************************/ -double GModelSpatialPointSource::eval(const GPhoton& photon) const +double GModelSpatialPointSource::eval(const GPhoton& photon, + const bool& gradients) const { // Set value dependent on source distance double value = (photon.dir().dist_deg(dir()) < tolerance) ? 1.0 : 0.0; @@ -296,27 +302,6 @@ double GModelSpatialPointSource::eval(const GPhoton& photon) const /***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] photon Incident photon. - * @return Model value. - * - * Evaluates the spatial part for a point source model and the gradient. - * It implements a delta function with respect to the coordinates of the - * source. For numerical reasons a certain tolerance is accepted (typically - * 0.1 arcsec, i.e. well below the angular resolution of gamma-ray - * telescopes). - * - * This method does not provide valid parameter gradients. - ***************************************************************************/ -double GModelSpatialPointSource::eval_gradients(const GPhoton& photon) const -{ - // Return value - return (eval(photon)); -} - - -/***********************************************************************//** * @brief Returns MC sky direction * * @param[in] energy Photon energy. diff --git a/src/model/GModelSpatialRadial.cpp b/src/model/GModelSpatialRadial.cpp index 97f9a23..7610776 100644 --- a/src/model/GModelSpatialRadial.cpp +++ b/src/model/GModelSpatialRadial.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialRadial.cpp - Abstract radial spatial model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2013 by Juergen Knoedlseder * + * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -162,38 +162,23 @@ GModelSpatialRadial& GModelSpatialRadial::operator=(const GModelSpatialRadial& m * @brief Return model value * * @param[in] photon Incident Photon. + * @param[in] gradients Compute gradients? + * @return Value of spatial radial model. * * Evaluates the radial spatial model value for a specific incident * @p photon. - ***************************************************************************/ -double GModelSpatialRadial::eval(const GPhoton& photon) const -{ - // Compute distance from source (in radians) - double theta = photon.dir().dist(dir()); - - // Evaluate model - double value = eval(theta, photon.energy(), photon.time()); - - // Return result - return value; -} - - -/***********************************************************************//** - * @brief Return model value and set analytical gradients - * - * @param[in] photon Incident Photon. * - * Evaluates the radial spatial model value and analytical model parameter - * gradients for a specific incident @p photon. + * If the @p gradients flag is true the method will also compute the + * parameter gradients for all model parameters. ***************************************************************************/ -double GModelSpatialRadial::eval_gradients(const GPhoton& photon) const +double GModelSpatialRadial::eval(const GPhoton& photon, + const bool& gradients) const { // Compute distance from source (in radians) double theta = photon.dir().dist(dir()); - // Evaluate model and set gradients - double value = eval_gradients(theta, photon.energy(), photon.time()); + // Evaluate model + double value = eval(theta, photon.energy(), photon.time(), gradients); // Return result return value; diff --git a/src/model/GModelSpatialRadialDisk.cpp b/src/model/GModelSpatialRadialDisk.cpp index 5033a56..24e1615 100644 --- a/src/model/GModelSpatialRadialDisk.cpp +++ b/src/model/GModelSpatialRadialDisk.cpp @@ -259,6 +259,7 @@ GModelSpatialRadialDisk* GModelSpatialRadialDisk::clone(void) const * @param[in] theta Angular distance from disk centre (radians). * @param[in] energy Photon energy. * @param[in] time Photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value. * * Evaluates the spatial component for a disk source model using @@ -281,10 +282,15 @@ GModelSpatialRadialDisk* GModelSpatialRadialDisk::clone(void) const * sky direction of interest, and * - \f${\tt m\_norm} = \frac{1}{2 \pi (1 - \cos r)} \f$ is a normalization * constant (see the update() method for details). + * + * The method will not compute analytical parameter gradients, even if the + * @p gradients argument is set to true. Radial disk parameter gradients + * need to be computed numerically. ***************************************************************************/ double GModelSpatialRadialDisk::eval(const double& theta, const GEnergy& energy, - const GTime& time) const + const GTime& time, + const bool& gradients) const { // Update precomputation cache update(); @@ -310,29 +316,6 @@ double GModelSpatialRadialDisk::eval(const double& theta, /***********************************************************************//** - * @brief Evaluate function and gradients (in units of sr^-1) - * - * @param[in] theta Angular distance from disk centre (radians). - * @param[in] energy Photon energy. - * @param[in] time Photon arrival time. - * @return Model value. - * - * Evaluates the function value. No gradient computation is implemented as - * radial models will be convolved with the instrument response and thus - * require the numerical computation of the derivatives. - * - * See the eval() method for more information. - ***************************************************************************/ -double GModelSpatialRadialDisk::eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const -{ - // Return value - return (eval(theta, energy, time)); -} - - -/***********************************************************************//** * @brief Return MC sky direction * * @param[in] energy Photon energy. diff --git a/src/model/GModelSpatialRadialGauss.cpp b/src/model/GModelSpatialRadialGauss.cpp index e6e3a35..1d43ff1 100644 --- a/src/model/GModelSpatialRadialGauss.cpp +++ b/src/model/GModelSpatialRadialGauss.cpp @@ -257,6 +257,7 @@ GModelSpatialRadialGauss* GModelSpatialRadialGauss::clone(void) const * @param[in] theta Angular distance from Gaussian centre (radians). * @param[in] energy Photon energy. * @param[in] time Photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value. * * Evaluates the spatial component for a Gaussian source model using @@ -271,12 +272,17 @@ GModelSpatialRadialGauss* GModelSpatialRadialGauss::clone(void) const * - \f$\theta\f$ is the angular separation from the centre of the model, and * - \f$\sigma\f$ is the Gaussian width. * + * The method will not compute analytical parameter gradients, even if the + * @p gradients argument is set to true. Radial disk parameter gradients + * need to be computed numerically. + * * @todo The Gaussian function is only correct in the small angle * approximation. ***************************************************************************/ double GModelSpatialRadialGauss::eval(const double& theta, const GEnergy& energy, - const GTime& time) const + const GTime& time, + const bool& gradients) const { // Compute value double sigma_rad = sigma() * gammalib::deg2rad; @@ -304,26 +310,6 @@ double GModelSpatialRadialGauss::eval(const double& theta, /***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] theta Angular distance from Gaussian centre (radians). - * @param[in] energy Photon energy. - * @param[in] time Photon arrival time. - * @return Model value. - * - * This method simply calls the eval() method as no analytical gradients - * will be computed. See the eval() method for more details. - ***************************************************************************/ -double GModelSpatialRadialGauss::eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const -{ - // Return value - return (eval(theta, energy, time)); -} - - -/***********************************************************************//** * @brief Returns MC sky direction * * @param[in] energy Photon energy. diff --git a/src/model/GModelSpatialRadialShell.cpp b/src/model/GModelSpatialRadialShell.cpp index 5fee237..5bf974a 100644 --- a/src/model/GModelSpatialRadialShell.cpp +++ b/src/model/GModelSpatialRadialShell.cpp @@ -258,11 +258,12 @@ GModelSpatialRadialShell* GModelSpatialRadialShell::clone(void) const /***********************************************************************//** - * @brief Evaluate function + * @brief Evaluate function (in units of sr^-1) * * @param[in] theta Angular distance from shell centre (radians). * @param[in] energy Photon energy. * @param[in] time Photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value. * * Evaluates the spatial part for a shell source model. The shell source @@ -290,10 +291,15 @@ GModelSpatialRadialShell* GModelSpatialRadialShell::clone(void) const * * Here, \f$\theta_{\rm in}\f$ and \f$\theta_{\rm out}\f$ are the shell inner * and outer radius. + * + * The method will not compute analytical parameter gradients, even if the + * @p gradients argument is set to true. Radial disk parameter gradients + * need to be computed numerically. ***************************************************************************/ double GModelSpatialRadialShell::eval(const double& theta, const GEnergy& energy, - const GTime& time) const + const GTime& time, + const bool& gradients) const { // Update precomputation cache update(); @@ -332,26 +338,6 @@ double GModelSpatialRadialShell::eval(const double& theta, /***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] theta Angular distance from shell centre (radians). - * @param[in] energy Photon energy. - * @param[in] time Photon arrival time. - * @return Model value. - * - * This method simply calls the eval() method as no analytical gradients will - * be computed. See the eval() method for details. - ***************************************************************************/ -double GModelSpatialRadialShell::eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const -{ - // Return value - return (eval(theta, energy, time)); -} - - -/***********************************************************************//** * @brief Returns MC sky direction * * @param[in] energy Photon energy. diff --git a/src/model/GModelSpectralBrokenPlaw.cpp b/src/model/GModelSpectralBrokenPlaw.cpp index b50e6e5..7e2856e 100644 --- a/src/model/GModelSpectralBrokenPlaw.cpp +++ b/src/model/GModelSpectralBrokenPlaw.cpp @@ -288,6 +288,7 @@ GModelSpectralBrokenPlaw* GModelSpectralBrokenPlaw::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value (ph/cm2/s/MeV). * * Evaluates @@ -302,60 +303,8 @@ GModelSpectralBrokenPlaw* GModelSpectralBrokenPlaw::clone(void) const * - \f${\tt m\_index}\f$ is the spectral index, and * - \f${\tt m\_breakenergy}\f$ is the breakenergy energy. * - * @todo The method expects that energy!=0. Otherwise Inf or NaN may result. - ***************************************************************************/ -double GModelSpectralBrokenPlaw::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Update the evaluation cache - update_eval_cache(srcEng); - - // Compute function value - double value = m_norm.value() * m_last_power; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralBrokenPlaw::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", m_norm=" << m_norm.value(); - std::cout << ", m_index1=" << m_index1.value(); - std::cout << ", m_breakenergy=" << m_breakenergy.value(); - std::cout << ", m_index2=" << m_index2.value(); - std::cout << ", m_last_power=" << m_last_power; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value (ph/cm2/s/MeV). - * - * Evaluates - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \left( \frac{E}{\tt m\_breakenergy} \right)^{\tt m\_index} - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization or prefactor, - * - \f${\tt m\_index}\f$ is the spectral index, and - * - \f${\tt m\_breakenergy}\f$ is the breakenergy energy. - * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = @@ -375,8 +324,9 @@ double GModelSpectralBrokenPlaw::eval(const GEnergy& srcEng, * * @todo The method expects that energy!=0. Otherwise Inf or NaN may result. ***************************************************************************/ -double GModelSpectralBrokenPlaw::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralBrokenPlaw::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { // Update the evaluation cache update_eval_cache(srcEng); @@ -384,46 +334,52 @@ double GModelSpectralBrokenPlaw::eval_gradients(const GEnergy& srcEng, // Compute function value double value = m_norm.value() * m_last_power; - // Compute normalisation gradient - double g_norm = (m_norm.is_free()) - ? m_norm.scale() * m_last_power : 0.0; - m_norm.factor_gradient(g_norm); - - // Compute index and break value gradients - if (srcEng.MeV() < m_breakenergy.value()) { - double g_index = (m_index1.is_free()) - ? value * m_index1.scale() * m_last_log_e_norm - : 0.0; - double g_break = (m_breakenergy.is_free()) - ? -value * m_last_index1 / m_breakenergy.factor_value() - : 0.0; - m_index1.factor_gradient(g_index); - m_index2.factor_gradient(0.0); - m_breakenergy.factor_gradient(g_break); - } - else { - double g_index = (m_index2.is_free()) - ? value * m_index2.scale() * m_last_log_e_norm - : 0.0; - double g_break = (m_breakenergy.is_free()) - ? -value * m_last_index2 / m_breakenergy.factor_value() - : 0.0; - m_index1.factor_gradient(0.0); - m_index2.factor_gradient(g_index); - m_breakenergy.factor_gradient(g_break); - } + // Optionally compute gradients + if (gradients) { + + // Compute normalisation gradient + double g_norm = (m_norm.is_free()) ? m_norm.scale() * m_last_power : 0.0; + m_norm.factor_gradient(g_norm); + + // Compute index and break value gradients + if (srcEng.MeV() < m_breakenergy.value()) { + double g_index = (m_index1.is_free()) + ? value * m_index1.scale() * m_last_log_e_norm + : 0.0; + double g_break = (m_breakenergy.is_free()) + ? -value * m_last_index1 / m_breakenergy.factor_value() + : 0.0; + m_index1.factor_gradient(g_index); + m_index2.factor_gradient(0.0); + m_breakenergy.factor_gradient(g_break); + } + else { + double g_index = (m_index2.is_free()) + ? value * m_index2.scale() * m_last_log_e_norm + : 0.0; + double g_break = (m_breakenergy.is_free()) + ? -value * m_last_index2 / m_breakenergy.factor_value() + : 0.0; + m_index1.factor_gradient(0.0); + m_index2.factor_gradient(g_index); + m_breakenergy.factor_gradient(g_break); + } + + } // 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)) { - std::cout << "*** ERROR: GModelSpectralBrokenPlaw::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralBrokenPlaw::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; - std::cout << ", e_norm=" << m_last_e_norm; - std::cout << ", power=" << m_last_power; - std::cout << ", g_norm=" << g_norm; + std::cout << ", m_norm=" << m_norm.value(); + std::cout << ", m_index1=" << m_index1.value(); + std::cout << ", m_breakenergy=" << m_breakenergy.value(); + std::cout << ", m_index2=" << m_index2.value(); + std::cout << ", m_last_power=" << m_last_power; std::cout << ")" << std::endl; } #endif @@ -872,7 +828,7 @@ void GModelSpectralBrokenPlaw::free_members(void) * * @param[in] energy Energy. * - * Updates the precomputation cache for eval() and eval_gradients() methods. + * Updates the precomputation cache for eval() the method. ***************************************************************************/ void GModelSpectralBrokenPlaw::update_eval_cache(const GEnergy& energy) const { diff --git a/src/model/GModelSpectralConst.cpp b/src/model/GModelSpectralConst.cpp index dfe2e33..06b7682 100644 --- a/src/model/GModelSpectralConst.cpp +++ b/src/model/GModelSpectralConst.cpp @@ -254,6 +254,7 @@ GModelSpectralConst* GModelSpectralConst::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value (ph/cm2/s/MeV). * * Evaluates @@ -265,54 +266,31 @@ GModelSpectralConst* GModelSpectralConst::clone(void) const * where * \f${\tt m\_norm}\f$ is the normalization constant in units of * ph/cm2/s/MeV. - ***************************************************************************/ -double GModelSpectralConst::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Compute function value - double value = m_norm.value(); - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate model value and gradient * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value (ph/cm2/s/MeV). - * - * Evaluates - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \f] - * - * where - * \f${\tt m\_norm}\f$ is the normalization constant in units of - * ph/cm2/s/MeV. - * - * The method also evaluates the partial derivative of the model with respect - * to the m_norm parameter using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = 1 * \f] ***************************************************************************/ -double GModelSpectralConst::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralConst::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { // Compute function value double value = m_norm.value(); - // Compute partial derivatives of the parameter values - double g_norm = (m_norm.is_free()) ? m_norm.scale() : 0.0; + // Optionally compute gradients + if (gradients) { + + // Compute partial derivatives of the parameter values + double g_norm = (m_norm.is_free()) ? m_norm.scale() : 0.0; + + // Set factor gradient + m_norm.factor_gradient(g_norm); - // Set factor gradient (the parameter gradient is obtained by dividing - // the factor gradient by the scale factor) - m_norm.factor_gradient(g_norm); + } // endif: gradient computation was requested // Return return value; diff --git a/src/model/GModelSpectralExpInvPlaw.cpp b/src/model/GModelSpectralExpInvPlaw.cpp index e952d5f..c59efcb 100644 --- a/src/model/GModelSpectralExpInvPlaw.cpp +++ b/src/model/GModelSpectralExpInvPlaw.cpp @@ -317,7 +317,8 @@ GModelSpectralExpInvPlaw* GModelSpectralExpInvPlaw::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. - * @return Model value. + * @param[in] gradients Compute gradients? + * @return Model value (ph/cm2/s/MeV). * * Evaluates * @@ -333,60 +334,8 @@ GModelSpectralExpInvPlaw* GModelSpectralExpInvPlaw::clone(void) const * - \f${\tt m\_index}\f$ is the spectral index, and * - \f${\tt m\_lambda}\f$ is the cut-off parameter. * - * @todo The method expects that pivot!=0. Otherwise Inf or NaN - * may result. We should add a test that prevents using invalid - * values. - ***************************************************************************/ -double GModelSpectralExpInvPlaw::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Update the evaluation cache - update_eval_cache(srcEng); - - // Compute function value - double value = m_norm.value() * m_last_power; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralExpPlaw::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", power=" << m_last_power; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value. - * - * Evaluates - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index} - * \exp \left( -{\tt m\_lambda}*E \right) - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization or prefactor, - * - \f${\tt m\_pivot}\f$ is the pivot energy, - * - \f${\tt m\_index}\f$ is the spectral index, and - * - \f${\tt m\_lambda}\f$ is the cut-off parameter. - * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = @@ -413,8 +362,9 @@ double GModelSpectralExpInvPlaw::eval(const GEnergy& srcEng, * may result. We should add a test that prevents using invalid * values. ***************************************************************************/ -double GModelSpectralExpInvPlaw::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralExpInvPlaw::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { // Update the evaluation cache update_eval_cache(srcEng); @@ -422,34 +372,41 @@ double GModelSpectralExpInvPlaw::eval_gradients(const GEnergy& srcEng, // Compute function value double value = m_norm.value() * m_last_power; - // Compute partial derivatives with respect to the parameter factor - // values. The partial derivatives with respect to the parameter - // values are obtained by division by the scale factor. - double g_norm = (m_norm.is_free()) - ? m_norm.scale() * m_last_power : 0.0; - double g_index = (m_index.is_free()) - ? value * m_index.scale() * std::log(m_last_e_norm) : 0.0; - double g_lambda = (m_lambda.is_free()) - ? -value * m_lambda.scale() * m_last_energy.MeV() : 0.0; - double g_pivot = (m_pivot.is_free()) - ? -value * m_last_index / m_pivot.factor_value() : 0.0; - - // Set gradients - m_norm.factor_gradient(g_norm); - m_index.factor_gradient(g_index); - m_lambda.factor_gradient(g_lambda); - m_pivot.factor_gradient(g_pivot); + // Optionally compute gradients + if (gradients) { + + // Compute partial derivatives with respect to the parameter factor + // values. The partial derivatives with respect to the parameter + // values are obtained by division by the scale factor. + double g_norm = (m_norm.is_free()) + ? m_norm.scale() * m_last_power + : 0.0; + double g_index = (m_index.is_free()) + ? value * m_index.scale() * std::log(m_last_e_norm) + : 0.0; + double g_lambda = (m_lambda.is_free()) + ? -value * m_lambda.scale() * m_last_energy.MeV() + : 0.0; + double g_pivot = (m_pivot.is_free()) + ? -value * m_last_index / m_pivot.factor_value() + : 0.0; + + // Set gradients + m_norm.factor_gradient(g_norm); + m_index.factor_gradient(g_index); + m_lambda.factor_gradient(g_lambda); + m_pivot.factor_gradient(g_pivot); + + } // 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)) { - std::cout << "*** ERROR: GModelSpectralExpPlaw::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralExpPlaw::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; - std::cout << ", e_norm=" << m_last_e_norm; - std::cout << ", e_lambda=" << m_last_e_lambda; std::cout << ", power=" << m_last_power; std::cout << ")" << std::endl; } @@ -875,7 +832,7 @@ void GModelSpectralExpInvPlaw::free_members(void) * * @param[in] energy Energy. * - * Updates the precomputation cache for eval() and eval_gradients() methods. + * Updates the precomputation cache for eval() method. ***************************************************************************/ void GModelSpectralExpInvPlaw::update_eval_cache(const GEnergy& energy) const { diff --git a/src/model/GModelSpectralExpPlaw.cpp b/src/model/GModelSpectralExpPlaw.cpp index 2d5eb46..5d548b0 100644 --- a/src/model/GModelSpectralExpPlaw.cpp +++ b/src/model/GModelSpectralExpPlaw.cpp @@ -287,7 +287,8 @@ GModelSpectralExpPlaw* GModelSpectralExpPlaw::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. - * @return Model value. + * @param[in] gradients Compute gradients? + * @return Model value (ph/cm2/s/MeV). * * Evaluates * @@ -303,60 +304,8 @@ GModelSpectralExpPlaw* GModelSpectralExpPlaw::clone(void) const * - \f${\tt m\_index}\f$ is the spectral index, and * - \f${\tt m\_ecut}\f$ is the cut off energy. * - * @todo The method expects that pivot!=0 and ecut!=0. Otherwise Inf or NaN - * may result. We should add a test that prevents using invalid - * values. - ***************************************************************************/ -double GModelSpectralExpPlaw::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Update the evaluation cache - update_eval_cache(srcEng); - - // Compute function value - double value = m_norm.value() * m_last_power; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralExpPlaw::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", power=" << m_last_power; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value. - * - * Evaluates - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index} - * \exp \left( \frac{-E}{\tt m\_ecut} \right) - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization or prefactor, - * - \f${\tt m\_index}\f$ is the spectral index, - * - \f${\tt m\_ecut}\f$ is the cut off energy, and - * - \f${\tt m\_pivot}\f$ is the pivot energy. - * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = @@ -383,8 +332,9 @@ double GModelSpectralExpPlaw::eval(const GEnergy& srcEng, * may result. We should add a test that prevents using invalid * values. ***************************************************************************/ -double GModelSpectralExpPlaw::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralExpPlaw::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { // Update the evaluation cache update_eval_cache(srcEng); @@ -392,34 +342,41 @@ double GModelSpectralExpPlaw::eval_gradients(const GEnergy& srcEng, // Compute function value double value = m_norm.value() * m_last_power; - // Compute partial derivatives with respect to the parameter factor - // values. The partial derivatives with respect to the parameter - // values are obtained by division by the scale factor. - double g_norm = (m_norm.is_free()) - ? m_norm.scale() * m_last_power : 0.0; - double g_index = (m_index.is_free()) - ? value * m_index.scale() * std::log(m_last_e_norm) : 0.0; - double g_ecut = (m_ecut.is_free()) - ? value * m_last_e_cut / m_ecut.factor_value() : 0.0; - double g_pivot = (m_pivot.is_free()) - ? -value * m_last_index / m_pivot.factor_value() : 0.0; - - // Set gradients - m_norm.factor_gradient(g_norm); - m_index.factor_gradient(g_index); - m_ecut.factor_gradient(g_ecut); - m_pivot.factor_gradient(g_pivot); + // Optionally compute gradients + if (gradients) { + + // Compute partial derivatives with respect to the parameter factor + // values. The partial derivatives with respect to the parameter + // values are obtained by division by the scale factor. + double g_norm = (m_norm.is_free()) + ? m_norm.scale() * m_last_power + : 0.0; + double g_index = (m_index.is_free()) + ? value * m_index.scale() * std::log(m_last_e_norm) + : 0.0; + double g_ecut = (m_ecut.is_free()) + ? value * m_last_e_cut / m_ecut.factor_value() + : 0.0; + double g_pivot = (m_pivot.is_free()) + ? -value * m_last_index / m_pivot.factor_value() + : 0.0; + + // Set gradients + m_norm.factor_gradient(g_norm); + m_index.factor_gradient(g_index); + m_ecut.factor_gradient(g_ecut); + m_pivot.factor_gradient(g_pivot); + + } // 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)) { - std::cout << "*** ERROR: GModelSpectralExpPlaw::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralExpPlaw::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; - std::cout << ", e_norm=" << m_last_e_norm; - std::cout << ", e_cut=" << m_last_e_cut; std::cout << ", power=" << m_last_power; std::cout << ")" << std::endl; } @@ -844,7 +801,7 @@ void GModelSpectralExpPlaw::free_members(void) * * @param[in] energy Energy. * - * Updates the precomputation cache for eval() and eval_gradients() methods. + * Updates the precomputation cache for eval() method. ***************************************************************************/ void GModelSpectralExpPlaw::update_eval_cache(const GEnergy& energy) const { diff --git a/src/model/GModelSpectralFunc.cpp b/src/model/GModelSpectralFunc.cpp index 6ec24f4..e837d34 100644 --- a/src/model/GModelSpectralFunc.cpp +++ b/src/model/GModelSpectralFunc.cpp @@ -234,52 +234,7 @@ GModelSpectralFunc* GModelSpectralFunc::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. - * @return Model value (ph/cm2/s/MeV). - * - * Evaluates - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization factor. - ***************************************************************************/ -double GModelSpectralFunc::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Interpolate function. This is done in log10-log10 space, but the - // linear value is returned. - double arg = m_log_nodes.interpolate(srcEng.log10MeV(), m_log_values); - double func = std::pow(10.0, arg); - - // Compute function value - double value = m_norm.value() * func; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralFunc::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", norm=" << norm(); - std::cout << ", func=" << func; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value (ph/cm2/s/MeV). * * Evaluates @@ -291,16 +246,17 @@ double GModelSpectralFunc::eval(const GEnergy& srcEng, * where * - \f${\tt m\_norm}\f$ is the normalization factor. * - * The method also evaluates the partial derivatives of the model with - * respect to the normalization parameter using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}} * \f] ***************************************************************************/ -double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralFunc::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { // Interpolate function. This is done in log10-log10 space, but the // linear value is returned. @@ -310,23 +266,27 @@ double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng, // Compute function value double value = m_norm.value() * func; - // Compute partial derivatives of the parameter values - double g_norm = (m_norm.is_free()) ? m_norm.scale() * func : 0.0; + // Optionally compute gradients + if (gradients) { - // Set gradients - m_norm.factor_gradient(g_norm); + // Compute partial derivatives of the parameter values + double g_norm = (m_norm.is_free()) ? m_norm.scale() * func : 0.0; + + // Set gradients + m_norm.factor_gradient(g_norm); + + } // 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)) { - std::cout << "*** ERROR: GModelSpectralFunc::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralFunc::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; std::cout << ", norm=" << norm(); std::cout << ", func=" << func; - std::cout << ", g_norm=" << g_norm; std::cout << ")" << std::endl; } #endif diff --git a/src/model/GModelSpectralGauss.cpp b/src/model/GModelSpectralGauss.cpp index 71ebaf0..fbb4fc3 100644 --- a/src/model/GModelSpectralGauss.cpp +++ b/src/model/GModelSpectralGauss.cpp @@ -236,57 +236,7 @@ GModelSpectralGauss* GModelSpectralGauss::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. - * @return Model value (ph/cm2/s/MeV). - * - * Evaluates - * - * \f[ - * S_{\rm E}(E | t) = \frac{m\_norm}{\sqrt{2\pi}m\_sigma} - * \exp(\frac{-(E-m\_mean)^2}{2 m\_sigma^2}) - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization, - * - \f${\tt m\_mean}\f$ is the mean energy, and - * - \f${\tt m\_sigma}\f$ is the energy width. - ***************************************************************************/ -double GModelSpectralGauss::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Get parameter values - double energy = srcEng.MeV(); - double norm = m_norm.value(); - double mean = m_mean.value(); - double sigma = m_sigma.value(); - - // Compute function value - double delta = energy - mean; - double term1 = (norm / sigma) * gammalib::inv_sqrt2pi; - double term2 = delta * delta / (2.0 * sigma * sigma); - double value = term1 * std::exp(-term2); - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralGauss::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate model value and gradient - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value (ph/cm2/s/MeV). * * Evaluates @@ -301,8 +251,8 @@ double GModelSpectralGauss::eval(const GEnergy& srcEng, * - \f${\tt m\_mean}\f$ is the mean energy, and * - \f${\tt m\_sigma}\f$ is the energy width. * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = @@ -319,10 +269,10 @@ double GModelSpectralGauss::eval(const GEnergy& srcEng, * \frac{S_{\rm E}(E | t)}{{\tt m\_sigma}} * \left( \frac{(E-m\_mean)^2}{m\_sigma^2} - 1 \right) * \f] - * ***************************************************************************/ -double GModelSpectralGauss::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralGauss::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { // Get parameter values double energy = srcEng.MeV(); @@ -343,25 +293,30 @@ double GModelSpectralGauss::eval_gradients(const GEnergy& srcEng, // Compute function value double value = term1 * eterm3; - // Compute partial derivatives with respect to the parameter factor - // values (partial differentials were determined analytically). - double g_norm = (m_norm.is_free()) - ? term2 * eterm3 * m_norm.scale() : 0.0; - double g_mean = (m_mean.is_free()) - ? value * term4 * m_mean.scale() : 0.0; - double g_sigma = (m_sigma.is_free()) - ? -term5 * eterm3 * (1.0 - (2.0 * term3)) * m_sigma.scale() - : 0.0; - - // Set gradients - m_norm.factor_gradient(g_norm); - m_mean.factor_gradient(g_mean); - m_sigma.factor_gradient(g_sigma); + // Optionally compute gradients + if (gradients) { + + // Compute partial derivatives with respect to the parameter factor + // values (partial differentials were determined analytically). + double g_norm = (m_norm.is_free()) + ? term2 * eterm3 * m_norm.scale() : 0.0; + double g_mean = (m_mean.is_free()) + ? value * term4 * m_mean.scale() : 0.0; + double g_sigma = (m_sigma.is_free()) + ? -term5 * eterm3 * (1.0 - (2.0 * term3)) * m_sigma.scale() + : 0.0; + + // Set gradients + m_norm.factor_gradient(g_norm); + m_mean.factor_gradient(g_mean); + m_sigma.factor_gradient(g_sigma); + + } // 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)) { - std::cout << "*** ERROR: GModelSpectralGauss::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralGauss::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; diff --git a/src/model/GModelSpectralLogParabola.cpp b/src/model/GModelSpectralLogParabola.cpp index 8906a9c..23e5e9e 100644 --- a/src/model/GModelSpectralLogParabola.cpp +++ b/src/model/GModelSpectralLogParabola.cpp @@ -296,6 +296,7 @@ GModelSpectralLogParabola* GModelSpectralLogParabola::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value (ph/cm2/s/MeV). * * Computes @@ -312,58 +313,8 @@ GModelSpectralLogParabola* GModelSpectralLogParabola::clone(void) const * - \f${\tt m\_curvature}\f$ is the spectral curvature, and * - \f${\tt m\_pivot}\f$ is the pivot energy. * - * @todo The method expects that energy!=0. Otherwise Inf or NaN may result. - ***************************************************************************/ -double GModelSpectralLogParabola::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Update the evaluation cache - update_eval_cache(srcEng); - - // Compute function value - double value = m_norm.value() * m_last_power; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralLogParabola::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", power=" << m_last_power; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value (ph/cm2/s/MeV). - * - * Computes - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \left( \frac{E}{\tt m\_pivot} \right)^{{\tt m\_index} + - * {\tt m\_curvature} \, \ln \frac{E}{\tt m\_pivot}} - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization or prefactor, - * - \f${\tt m\_index}\f$ is the spectral index, - * - \f${\tt m\_curvature}\f$ is the spectral curvature, and - * - \f${\tt m\_pivot}\f$ is the pivot energy. - * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = @@ -389,8 +340,9 @@ double GModelSpectralLogParabola::eval(const GEnergy& srcEng, * * @todo The method expects that energy!=0. Otherwise Inf or NaN may result. ***************************************************************************/ -double GModelSpectralLogParabola::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralLogParabola::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { // Update the evaluation cache update_eval_cache(srcEng); @@ -398,38 +350,38 @@ double GModelSpectralLogParabola::eval_gradients(const GEnergy& srcEng, // Compute function value double value = m_norm.value() * m_last_power; - // Compute partial derivatives of the parameter values - double g_norm = (m_norm.is_free()) - ? m_norm.scale() * m_last_power : 0.0; - double g_index = (m_index.is_free()) - ? value * m_index.scale() * m_last_log_e_norm : 0.0; - double g_curvature = (m_curvature.is_free()) - ? value * m_curvature.scale() * m_last_log_e_norm * - m_last_log_e_norm : 0.0; - double g_pivot = (m_pivot.is_free()) - ? -value * (m_last_exponent + m_curvature.value() * - m_last_log_e_norm) / m_pivot.factor_value() : 0.0; - - // Set gradients - m_norm.factor_gradient(g_norm); - m_index.factor_gradient(g_index); - m_curvature.factor_gradient(g_curvature); - m_pivot.factor_gradient(g_pivot); + // Optionally compute gradients + if (gradients) { + + // Compute partial derivatives of the parameter values + double g_norm = (m_norm.is_free()) + ? m_norm.scale() * m_last_power : 0.0; + double g_index = (m_index.is_free()) + ? value * m_index.scale() * m_last_log_e_norm : 0.0; + double g_curvature = (m_curvature.is_free()) + ? value * m_curvature.scale() * m_last_log_e_norm * + m_last_log_e_norm : 0.0; + double g_pivot = (m_pivot.is_free()) + ? -value * (m_last_exponent + m_curvature.value() * + m_last_log_e_norm) / m_pivot.factor_value() : 0.0; + + // Set gradients + m_norm.factor_gradient(g_norm); + m_index.factor_gradient(g_index); + m_curvature.factor_gradient(g_curvature); + m_pivot.factor_gradient(g_pivot); + + } // 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)) { - std::cout << "*** ERROR: GModelSpectralLogParabola::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralLogParabola::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; - std::cout << ", e_norm=" << m_last_e_norm; std::cout << ", power=" << m_last_power; - std::cout << ", g_norm=" << g_norm; - std::cout << ", g_index=" << g_index; - std::cout << ", g_curvature=" << g_curvature; - std::cout << ", g_pivot=" << g_pivot; std::cout << ")" << std::endl; } #endif @@ -884,7 +836,7 @@ void GModelSpectralLogParabola::free_members(void) * * @param[in] energy Energy. * - * Updates the precomputation cache for eval() and eval_gradients() methods. + * Updates the precomputation cache for eval() method. ***************************************************************************/ void GModelSpectralLogParabola::update_eval_cache(const GEnergy& energy) const { diff --git a/src/model/GModelSpectralNodes.cpp b/src/model/GModelSpectralNodes.cpp index ac216a4..1338c63 100644 --- a/src/model/GModelSpectralNodes.cpp +++ b/src/model/GModelSpectralNodes.cpp @@ -210,6 +210,7 @@ GModelSpectralNodes* GModelSpectralNodes::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value (ph/cm2/s/MeV). * * Computes @@ -233,75 +234,18 @@ GModelSpectralNodes* GModelSpectralNodes::clone(void) const * and * \f$(\log {\tt m\_energies[i+1]}, \log{\tt m\_values[i+1]})\f$ * to the requested energy \f$\log E\f$. - ***************************************************************************/ -double GModelSpectralNodes::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Update evaluation cache - update_eval_cache(); - - // Interpolate function. This is done in log10-log10 space, but the - // linear value is returned. - double exponent = m_log_energies.interpolate(srcEng.log10MeV(), m_log_values); - double value = std::pow(10.0, exponent); - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralNodes::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", exponent=" << exponent; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value (ph/cm2/s/MeV). * - * Computes - * - * \f[ - * S_{\rm E}(E | t) = - * 10^{(\log {\tt m\_values[i]}) w_{i} + - * (\log {\tt m\_values[i+1]}) w_{i+1}} - * \f] - * - * where - * - \f${\tt m\_values[i]}\f$ is the intensity of node \f$i\f$, - * - \f${\tt m\_values[i+1]}\f$ is the intensity of node \f$i+1\f$, - * - \f$w_{i}\f$ is the weighting of node \f$i\f$, - * - \f$w_{i+1}\f$ is the weighting of node \f$i+1\f$, and - * - \f${\tt m\_energies[i]} <= E <= {\tt m\_energies[i+1]}\f$. - * - * The weightings \f$w_{i}\f$ and \f$w_{i+1}\f$ are computed by linear - * interpolation (in the log-log plane) between the nodes - * \f$(\log {\tt m\_energies[i]}, \log{\tt m\_values[i]})\f$ - * and - * \f$(\log {\tt m\_energies[i+1]}, \log{\tt m\_values[i+1]})\f$ - * to the requested energy \f$\log E\f$. - * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using + * If the @p gradients flag is true the method will also compute the partial + * derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_values[i]}} = * \frac{S_{\rm E}(E | t) \, w_i}{{\tt m\_values[i]}} * \f] ***************************************************************************/ -double GModelSpectralNodes::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralNodes::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { // Update evaluation cache update_eval_cache(); @@ -322,27 +266,32 @@ double GModelSpectralNodes::eval_gradients(const GEnergy& srcEng, // Compute linear value double value = std::pow(10.0, exponent); - // Initialise gradients - for (int i = 0; i < m_values.size(); ++i) { - m_values[i].factor_gradient(0.0); - } + // Optionally compute partial derivatives + if (gradients) { - // Gradient for left node - if (m_values[inx_left].is_free()) { - double grad = value * wgt_left / m_values[inx_left].factor_value(); - m_values[inx_left].factor_gradient(grad); - } + // Initialise gradients + for (int i = 0; i < m_values.size(); ++i) { + m_values[i].factor_gradient(0.0); + } - // Gradient for right node - if (m_values[inx_right].is_free()) { - double grad = value * wgt_right / m_values[inx_right].factor_value(); - m_values[inx_right].factor_gradient(grad); - } + // Gradient for left node + if (m_values[inx_left].is_free()) { + double grad = value * wgt_left / m_values[inx_left].factor_value(); + m_values[inx_left].factor_gradient(grad); + } + + // Gradient for right node + if (m_values[inx_right].is_free()) { + double grad = value * wgt_right / m_values[inx_right].factor_value(); + m_values[inx_right].factor_gradient(grad); + } + + } // endif: computed partial derivatives // Compile option: Check for NaN/Inf #if defined(G_NAN_CHECK) if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralNodes::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralNodes::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; diff --git a/src/model/GModelSpectralPlaw.cpp b/src/model/GModelSpectralPlaw.cpp index 463c41d..feb1e78 100644 --- a/src/model/GModelSpectralPlaw.cpp +++ b/src/model/GModelSpectralPlaw.cpp @@ -276,6 +276,7 @@ GModelSpectralPlaw* GModelSpectralPlaw::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value (ph/cm2/s/MeV). * * Evaluates @@ -290,59 +291,8 @@ GModelSpectralPlaw* GModelSpectralPlaw::clone(void) const * - \f${\tt m\_index}\f$ is the spectral index, and * - \f${\tt m\_pivot}\f$ is the pivot energy. * - * @todo The method expects that energy!=0. Otherwise Inf or NaN may result. - ***************************************************************************/ -double GModelSpectralPlaw::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Update the evaluation cache - update_eval_cache(srcEng); - - // Compute function value - double value = m_norm.value() * m_last_power; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralPlaw::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", m_norm=" << m_norm.value(); - std::cout << ", m_index=" << m_index.value(); - std::cout << ", m_pivot=" << m_pivot.value(); - std::cout << ", m_last_power=" << m_last_power; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value (ph/cm2/s/MeV). - * - * Evaluates - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index} - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization or prefactor, - * - \f${\tt m\_index}\f$ is the spectral index, and - * - \f${\tt m\_pivot}\f$ is the pivot energy. - * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = @@ -362,8 +312,9 @@ double GModelSpectralPlaw::eval(const GEnergy& srcEng, * * @todo The method expects that energy!=0. Otherwise Inf or NaN may result. ***************************************************************************/ -double GModelSpectralPlaw::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralPlaw::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { // Update the evaluation cache update_eval_cache(srcEng); @@ -371,32 +322,36 @@ double GModelSpectralPlaw::eval_gradients(const GEnergy& srcEng, // Compute function value double value = m_norm.value() * m_last_power; - // Compute partial derivatives of the parameter values - double g_norm = (m_norm.is_free()) - ? m_norm.scale() * m_last_power : 0.0; - double g_index = (m_index.is_free()) - ? value * m_index.scale() * m_last_log_e_norm : 0.0; - double g_pivot = (m_pivot.is_free()) - ? -value * m_last_index / m_pivot.factor_value() : 0.0; + // Optionally compute gradients + if (gradients) { - // Set gradients - m_norm.factor_gradient(g_norm); - m_index.factor_gradient(g_index); - m_pivot.factor_gradient(g_pivot); + // Compute partial derivatives of the parameter values + double g_norm = (m_norm.is_free()) + ? m_norm.scale() * m_last_power : 0.0; + double g_index = (m_index.is_free()) + ? value * m_index.scale() * m_last_log_e_norm : 0.0; + double g_pivot = (m_pivot.is_free()) + ? -value * m_last_index / m_pivot.factor_value() : 0.0; + + // Set gradients + m_norm.factor_gradient(g_norm); + m_index.factor_gradient(g_index); + m_pivot.factor_gradient(g_pivot); + + } // 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)) { - std::cout << "*** ERROR: GModelSpectralPlaw::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralPlaw::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; - std::cout << ", e_norm=" << m_last_e_norm; - std::cout << ", power=" << m_last_power; - std::cout << ", g_norm=" << g_norm; - std::cout << ", g_index=" << g_index; - std::cout << ", g_pivot=" << g_pivot; + std::cout << ", m_norm=" << m_norm.value(); + std::cout << ", m_index=" << m_index.value(); + std::cout << ", m_pivot=" << m_pivot.value(); + std::cout << ", m_last_power=" << m_last_power; std::cout << ")" << std::endl; } #endif @@ -770,7 +725,7 @@ void GModelSpectralPlaw::free_members(void) * * @param[in] energy Energy. * - * Updates the precomputation cache for eval() and eval_gradients() methods. + * Updates the precomputation cache for eval() method. ***************************************************************************/ void GModelSpectralPlaw::update_eval_cache(const GEnergy& energy) const { diff --git a/src/model/GModelSpectralPlawEnergyFlux.cpp b/src/model/GModelSpectralPlawEnergyFlux.cpp index 5d73ad0..ce58aaa 100644 --- a/src/model/GModelSpectralPlawEnergyFlux.cpp +++ b/src/model/GModelSpectralPlawEnergyFlux.cpp @@ -279,6 +279,7 @@ GModelSpectralPlawEnergyFlux* GModelSpectralPlawEnergyFlux::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value (ph/cm2/s/MeV). * * Computes @@ -306,71 +307,9 @@ GModelSpectralPlawEnergyFlux* GModelSpectralPlawEnergyFlux::clone(void) const * - \f${\tt m\_eflux}\f$ is the energy flux between * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and * - \f${\tt m\_index}\f$ is the spectral index. - ***************************************************************************/ -double GModelSpectralPlawEnergyFlux::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Update precomputed values - update(srcEng); - - // Compute function value - double value = m_eflux.value() * gammalib::erg2MeV * m_norm * m_power; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralPlawEnergyFlux::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", eflux=" << eflux(); - std::cout << ", m_norm=" << m_norm; - std::cout << ", m_power=" << m_power; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value (ph/cm2/s/MeV). - * - * Computes - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_eflux} - * \frac{{\tt m\_index}+2} - * {{\tt e\_max}^{{\tt m\_index}+2} - - * {\tt e\_min}^{{\tt m\_index}+2}} - * E^{\tt m\_index} - * \f] - * - * for \f${\tt m\_index} \ne -2\f$ and - * - * \f[ - * S_{\rm E}(E | t) = - * \frac{{\tt m\_eflux}} - * {\log {\tt e\_max} - \log {\tt e\_min}} - * E^{\tt m\_index} - * \f] - * - * for \f${\tt m\_index} = -2\f$, where - * - \f${\tt e\_min}\f$ is the minimum energy of an interval, - * - \f${\tt e\_max}\f$ is the maximum energy of an interval, - * - \f${\tt m\_eflux}\f$ is the energy flux between - * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and - * - \f${\tt m\_index}\f$ is the spectral index. - * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_eflux}} = @@ -399,38 +338,38 @@ double GModelSpectralPlawEnergyFlux::eval(const GEnergy& srcEng, * * No partial derivatives are supported for the energy boundaries. ***************************************************************************/ -double GModelSpectralPlawEnergyFlux::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralPlawEnergyFlux::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { - // Initialise gradients - double g_eflux = 0.0; - double g_index = 0.0; - // Update precomputed values update(srcEng); // Compute function value double value = m_eflux.value() * gammalib::erg2MeV * m_norm * m_power; - // Integral flux gradient - if (m_eflux.is_free() && m_eflux.factor_value() > 0.0) { - g_eflux = value / m_eflux.factor_value(); - } + // Optionally compute gradients + if (gradients) { - // Index gradient - if (m_index.is_free()) { - g_index = value * (m_g_norm + gammalib::ln10 * srcEng.log10MeV()) * - m_index.scale(); - } + // Compute partial derivatives of the parameter values + double g_eflux = (m_eflux.is_free() && m_eflux.factor_value() > 0.0) + ? g_eflux = value / m_eflux.factor_value() + : 0.0; + double g_index = (m_index.is_free()) + ? value * (m_g_norm + gammalib::ln10 * + srcEng.log10MeV()) * m_index.scale() + : 0.0; - // Set gradients - m_eflux.factor_gradient(g_eflux); - m_index.factor_gradient(g_index); + // Set gradients + m_eflux.factor_gradient(g_eflux); + m_index.factor_gradient(g_index); + + } // 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)) { - std::cout << "*** ERROR: GModelSpectralPlawEnergyFlux::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralPlawEnergyFlux::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; @@ -438,8 +377,6 @@ double GModelSpectralPlawEnergyFlux::eval_gradients(const GEnergy& srcEng, std::cout << ", eflux=" << eflux(); std::cout << ", m_norm=" << m_norm; std::cout << ", m_power=" << m_power; - std::cout << ", g_eflux=" << g_eflux; - std::cout << ", g_index=" << g_index; std::cout << ")" << std::endl; } #endif diff --git a/src/model/GModelSpectralPlawPhotonFlux.cpp b/src/model/GModelSpectralPlawPhotonFlux.cpp index 3f6d72b..d032b8a 100644 --- a/src/model/GModelSpectralPlawPhotonFlux.cpp +++ b/src/model/GModelSpectralPlawPhotonFlux.cpp @@ -289,6 +289,7 @@ GModelSpectralPlawPhotonFlux* GModelSpectralPlawPhotonFlux::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value (ph/cm2/s/MeV). * * Computes @@ -316,71 +317,9 @@ GModelSpectralPlawPhotonFlux* GModelSpectralPlawPhotonFlux::clone(void) const * - \f${\tt m\_flux}\f$ is the photon flux between * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and * - \f${\tt m\_index}\f$ is the spectral index. - ***************************************************************************/ -double GModelSpectralPlawPhotonFlux::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Update precomputed values - update(srcEng); - - // Compute function value - double value = m_flux.value() * m_norm * m_power; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralPlawPhotonFlux::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", flux=" << flux(); - std::cout << ", m_norm=" << m_norm; - std::cout << ", m_power=" << m_power; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value (ph/cm2/s/MeV). * - * Computes - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_flux} - * \frac{{\tt m\_index}+1} - * {{\tt e\_max}^{{\tt m\_index}+1} - - * {\tt e\_min}^{{\tt m\_index}+1}} - * E^{\tt m\_index} - * \f] - * - * for \f${\tt m\_index} \ne -1\f$ and - * - * \f[ - * S_{\rm E}(E | t) = - * \frac{{\tt m\_flux}} - * {\log {\tt e\_max} - \log {\tt e\_min}} - * E^{\tt m\_index} - * \f] - * - * for \f${\tt m\_index} = -1\f$, where - * - \f${\tt e\_min}\f$ is the minimum energy of an interval, - * - \f${\tt e\_max}\f$ is the maximum energy of an interval, - * - \f${\tt m\_flux}\f$ is the photon flux between - * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and - * - \f${\tt m\_index}\f$ is the spectral index. - * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_flux}} = @@ -406,41 +345,38 @@ double GModelSpectralPlawPhotonFlux::eval(const GEnergy& srcEng, * \f] * * for \f${\tt m\_index} = -1\f$. - * - * No partial derivatives are supported for the energy boundaries. ***************************************************************************/ -double GModelSpectralPlawPhotonFlux::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralPlawPhotonFlux::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { - // Initialise gradients - double g_flux = 0.0; - double g_index = 0.0; - // Update precomputed values update(srcEng); // Compute function value double value = m_flux.value() * m_norm * m_power; - // Photon flux gradient - if (m_flux.is_free() && m_flux.factor_value() > 0.0) { - g_flux = value / m_flux.factor_value(); - } + // Optionally compute gradients + if (gradients) { - // Index gradient - if (m_index.is_free()) { - g_index = value * (m_g_norm + gammalib::ln10 * srcEng.log10MeV()) * - m_index.scale(); - } + // Compute partial derivatives of the parameter values + double g_flux = (m_flux.is_free() && m_flux.factor_value() > 0.0) + ? value / m_flux.factor_value() + : 0.0; + double g_index = (m_index.is_free()) + ? value * (m_g_norm + gammalib::ln10 * + srcEng.log10MeV()) * m_index.scale() : 0.0; - // Set gradients - m_flux.factor_gradient(g_flux); - m_index.factor_gradient(g_index); + // Set gradients + m_flux.factor_gradient(g_flux); + m_index.factor_gradient(g_index); + + } // 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)) { - std::cout << "*** ERROR: GModelSpectralPlawPhotonFlux::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralPlawPhotonFlux::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; @@ -448,8 +384,6 @@ double GModelSpectralPlawPhotonFlux::eval_gradients(const GEnergy& srcEng, std::cout << ", flux=" << flux(); std::cout << ", m_norm=" << m_norm; std::cout << ", m_power=" << m_power; - std::cout << ", g_flux=" << g_flux; - std::cout << ", g_index=" << g_index; std::cout << ")" << std::endl; } #endif diff --git a/src/model/GModelSpectralSuperExpPlaw.cpp b/src/model/GModelSpectralSuperExpPlaw.cpp index 2b23601..e9e3015 100644 --- a/src/model/GModelSpectralSuperExpPlaw.cpp +++ b/src/model/GModelSpectralSuperExpPlaw.cpp @@ -297,7 +297,8 @@ GModelSpectralSuperExpPlaw* GModelSpectralSuperExpPlaw::clone(void) const * * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. - * @return Model value. + * @param[in] gradients Compute gradients? + * @return Model value (ph/cm2/s/MeV). * * Evaluates * @@ -315,62 +316,8 @@ GModelSpectralSuperExpPlaw* GModelSpectralSuperExpPlaw::clone(void) const * - \f${\tt m\_ecut}\f$ is the cut off energy, and * - \f${\tt m\_index2}\f$ is the index of the cut off * - * @todo The method expects that pivot!=0 and ecut!=0. Otherwise Inf or NaN - * may result. We should add a test that prevents using invalid - * values. - ***************************************************************************/ -double GModelSpectralSuperExpPlaw::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Update the evaluation cache - update_eval_cache(srcEng); - - // Compute function value - double value = m_norm.value() * m_last_power; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralSuperExpPlaw::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", power=" << m_last_power; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value. - * - * Evaluates - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index1} - * \exp \left( - \left(\frac{-E}{\tt m\_ecut} \right)^{\tt m_index2} - * \right) - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization or prefactor, - * - \f${\tt m\_index1}\f$ is the spectral index, - * - \f${\tt m\_ecut}\f$ is the cut off energy - * - \f${\tt m\_pivot}\f$ is the pivot energy, and - * - \f${\tt m\_index2}\f$ is the index of the cut off - * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using + * If the @p gradients flag is true the method will also compute the + * partial derivatives of the model with respect to the parameters using * * \f[ * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = @@ -404,8 +351,9 @@ double GModelSpectralSuperExpPlaw::eval(const GEnergy& srcEng, * may result. We should add a test that prevents using invalid * values. ***************************************************************************/ -double GModelSpectralSuperExpPlaw::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) +double GModelSpectralSuperExpPlaw::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const { // Update the evaluation cache update_eval_cache(srcEng); @@ -413,42 +361,47 @@ double GModelSpectralSuperExpPlaw::eval_gradients(const GEnergy& srcEng, // Compute function value double value = m_norm.value() * m_last_power; - // Compute partial derivatives with respect to the parameter factor - // values. The partial derivatives with respect to the parameter - // values are obtained by multiplication with the scale factor. - double g_norm = (m_norm.is_free()) - ? m_norm.scale() * m_last_power - : 0.0; - double g_index1 = (m_index1.is_free()) - ? value * m_index1.scale() * std::log(m_last_e_norm) - : 0.0; - double g_ecut = (m_ecut.is_free()) - ? value * m_last_index2 * m_last_exponent / m_ecut.factor_value() - : 0.0; - double g_pivot = (m_pivot.is_free()) - ? -value * m_last_index1 / m_pivot.factor_value() - : 0.0; - double g_index2 = (m_index2.is_free()) - ? -value * m_index2.scale() * std::log(m_last_e_cut) * m_last_exponent - : 0.0; - - // Set gradients - m_norm.factor_gradient(g_norm); - m_index1.factor_gradient(g_index1); - m_ecut.factor_gradient(g_ecut); - m_pivot.factor_gradient(g_pivot); - m_index2.factor_gradient(g_index2); + // Optionally compute gradients + if (gradients) { + + // Compute partial derivatives with respect to the parameter factor + // values. The partial derivatives with respect to the parameter + // values are obtained by multiplication with the scale factor. + double g_norm = (m_norm.is_free()) + ? m_norm.scale() * m_last_power + : 0.0; + double g_index1 = (m_index1.is_free()) + ? value * m_index1.scale() * std::log(m_last_e_norm) + : 0.0; + double g_ecut = (m_ecut.is_free()) + ? value * m_last_index2 * m_last_exponent / + m_ecut.factor_value() + : 0.0; + double g_pivot = (m_pivot.is_free()) + ? -value * m_last_index1 / m_pivot.factor_value() + : 0.0; + double g_index2 = (m_index2.is_free()) + ? -value * m_index2.scale() * + std::log(m_last_e_cut) * m_last_exponent + : 0.0; + + // Set gradients + m_norm.factor_gradient(g_norm); + m_index1.factor_gradient(g_index1); + m_ecut.factor_gradient(g_ecut); + m_pivot.factor_gradient(g_pivot); + m_index2.factor_gradient(g_index2); + + } // 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)) { - std::cout << "*** ERROR: GModelSpectralSuperExpPlaw::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralSuperExpPlaw::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; - std::cout << ", e_norm=" << m_last_e_norm; - std::cout << ", e_cut=" << m_last_e_cut; std::cout << ", power=" << m_last_power; std::cout << ")" << std::endl; } @@ -894,7 +847,7 @@ void GModelSpectralSuperExpPlaw::free_members(void) * * @param[in] energy Energy. * - * Updates the precomputation cache for eval() and eval_gradients() methods. + * Updates the precomputation cache for eval() method. ***************************************************************************/ void GModelSpectralSuperExpPlaw::update_eval_cache(const GEnergy& energy) const { diff --git a/src/model/GModelTemporalConst.cpp b/src/model/GModelTemporalConst.cpp index 3e9bf05..8bccda5 100644 --- a/src/model/GModelTemporalConst.cpp +++ b/src/model/GModelTemporalConst.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelTemporalConst.cpp - Temporal constant model class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2015 by Juergen Knoedlseder * + * copyright (C) 2009-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -225,6 +225,8 @@ GModelTemporalConst* GModelTemporalConst::clone(void) const * @brief Evaluate function * * @param[in] srcTime True photon arrival time (not used). + * @param[in] gradients Compute gradients? + * @return Value to temporal model. * * Computes * @@ -234,49 +236,30 @@ GModelTemporalConst* GModelTemporalConst::clone(void) const * * where * \f${\tt m\_norm}\f$ is the normalization constant. - ***************************************************************************/ -double GModelTemporalConst::eval(const GTime& srcTime) const -{ - // Compute function value - double value = norm(); - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients * - * @param[in] srcTime True photon arrival time (not used). - * - * Computes - * - * \f[ - * S_{\rm t}(t) = {\tt m\_norm} - * \f] - * - * where - * \f${\tt m\_norm}\f$ is the normalization constant. - * - * The method also evaluates the partial derivatives of the model with - * respect to the normalization parameter using + * If the @p gradients flag is true the method will also evaluate the partial + * derivatives of the model with respect to the normalization parameter using * * \f[ * \frac{\delta S_{\rm t}(t)}{\delta {\tt m\_norm}} = 1 * \f] ***************************************************************************/ -double GModelTemporalConst::eval_gradients(const GTime& srcTime) +double GModelTemporalConst::eval(const GTime& srcTime, + const bool& gradients) const { // Compute function value double value = m_norm.value(); - // Compute partial derivatives of the parameter values - double g_norm = (m_norm.is_free()) ? m_norm.scale() : 0.0; + // Optionally compute partial derivatives + if (gradients) { + + // Compute partial derivatives of the parameter values + double g_norm = (m_norm.is_free()) ? m_norm.scale() : 0.0; + + // Set factor gradient + m_norm.factor_gradient(g_norm); - // Set factor gradient (the parameter gradient is obtained by dividing - // the factor gradient by the scale factor) - m_norm.factor_gradient(g_norm); + } // endif: computed partial derivatives // Return return value; diff --git a/src/model/GModels.cpp b/src/model/GModels.cpp index d0c4cd8..1cc6765 100644 --- a/src/model/GModels.cpp +++ b/src/model/GModels.cpp @@ -853,11 +853,18 @@ GOptimizerPars GModels::pars(void) const * * @param[in] event Observed event. * @param[in] obs Observation. + * @param[in] gradients Compute gradients? + * @return Value of models. * * Evaluates the sum of all models for the specified event and observation. * Only valid models are considered in this evaluation. + * + * If the @p gradients flag is true the method will also set the parameter + * gradients of the model parameters. ***************************************************************************/ -double GModels::eval(const GEvent& event, const GObservation& obs) const +double GModels::eval(const GEvent& event, + const GObservation& obs, + const bool& gradients) const { // Initialise function value double value = 0.0; @@ -865,7 +872,7 @@ double GModels::eval(const GEvent& event, const GObservation& obs) const // Evaluate function for all models for (int i = 0; i < size(); ++i) { if (m_models[i]->is_valid(obs.instrument(), obs.id())) { - value += m_models[i]->eval(event, obs); + value += m_models[i]->eval(event, obs, gradients); } } @@ -875,31 +882,6 @@ double GModels::eval(const GEvent& event, const GObservation& obs) const /***********************************************************************//** - * @brief Evaluate sum and gradients of all models - * - * @param[in] event Observed event. - * @param[in] obs Observation. - * - * Evaluates the sum and the parameter gradients of all models for the - * specified event and observation. - ***************************************************************************/ -double GModels::eval_gradients(const GEvent& event, - const GObservation& obs) const -{ - // Initialise function value - double value = 0.0; - - // Evaluate function for all models - for (int i = 0; i < size(); ++i) { - value += m_models[i]->eval_gradients(event, obs); - } - - // Return - return value; -} - - -/***********************************************************************//** * @brief Print models * * @param[in] chatter Chattiness (defaults to NORMAL). diff --git a/src/obs/GObservation.cpp b/src/obs/GObservation.cpp index bef3481..8ba4e48 100644 --- a/src/obs/GObservation.cpp +++ b/src/obs/GObservation.cpp @@ -290,12 +290,7 @@ double GObservation::model(const GModels& models, // is used, don't compute model gradients as we cannot // use them. This is somehow a kluge, but makes the // code faster - if (response()->use_edisp()) { - model += mptr->eval(event, *this); - } - else { - model += mptr->eval_gradients(event, *this); - } + model += mptr->eval(event, *this, !response()->use_edisp()); // Optionally determine model gradients. If the model has a // gradient then use it, unless we have energy dispersion. diff --git a/src/obs/GResponse.cpp b/src/obs/GResponse.cpp index 1d17755..74f55a7 100644 --- a/src/obs/GResponse.cpp +++ b/src/obs/GResponse.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GResponse.cpp - Abstract response base class * * ----------------------------------------------------------------------- * - * copyright (C) 2008-2015 by Juergen Knoedlseder * + * copyright (C) 2008-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -332,28 +332,17 @@ double GResponse::eval_prob(const GModelSky& model, irf *= model.scale(obs.instrument()).value(); } - // Case A: evaluate gradients - if (grad) { + // Evaluate spectral and temporal components + double spec = (model.spectral() != NULL) + ? model.spectral()->eval(srcEng, srcTime, grad) : 1.0; + double temp = (model.temporal() != NULL) + ? model.temporal()->eval(srcTime, grad) : 1.0; - // Evaluate source model - double spec = (model.spectral() != NULL) ? model.spectral()->eval_gradients(srcEng, srcTime) : 1.0; - double temp = (model.temporal() != NULL) ? model.temporal()->eval_gradients(srcTime) : 1.0; - - // Set probability - prob = spec * temp * irf; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(prob) || gammalib::is_infinite(prob)) { - std::cout << "*** ERROR: GResponse::eval_prob:"; - std::cout << " NaN/Inf encountered"; - std::cout << " (prob=" << prob; - std::cout << ", spec=" << spec; - std::cout << ", temp=" << temp; - std::cout << ", irf=" << irf; - std::cout << ")" << std::endl; - } - #endif + // Compute probability + prob = spec * temp * irf; + + // Optionally compute partial derivatives + if (grad) { // Multiply factors to spectral gradients if (model.spectral() != NULL) { @@ -375,32 +364,20 @@ double GResponse::eval_prob(const GModelSky& model, } } - } // endif: gradient evaluation has been requested - - // Case B: evaluate no gradients - else { - - // Evaluate source model - double spec = (model.spectral() != NULL) ? model.spectral()->eval(srcEng, srcTime) : 1.0; - double temp = (model.temporal() != NULL) ? model.temporal()->eval(srcTime) : 1.0; - - // Set probability - prob = spec * temp * irf; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(prob) || gammalib::is_infinite(prob)) { - std::cout << "*** ERROR: GResponse::eval_prob:"; - std::cout << " NaN/Inf encountered"; - std::cout << " (prob=" << prob; - std::cout << ", spec=" << spec; - std::cout << ", temp=" << temp; - std::cout << ", irf=" << irf; - std::cout << ")" << std::endl; - } - #endif + } // endif: computed partial derivatives + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(prob) || gammalib::is_infinite(prob)) { + std::cout << "*** ERROR: GResponse::eval_prob:"; + std::cout << " NaN/Inf encountered"; + std::cout << " (prob=" << prob; + std::cout << ", spec=" << spec; + std::cout << ", temp=" << temp; + std::cout << ", irf=" << irf; + std::cout << ")" << std::endl; } + #endif } // endif: Gamma-ray source model had a spatial component diff --git a/test/testinst/GTestModelData.hpp b/test/testinst/GTestModelData.hpp index 19db91e..860288a 100644 --- a/test/testinst/GTestModelData.hpp +++ b/test/testinst/GTestModelData.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GTestModelData.hpp - Test data model class * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2014 by Jean-Baptiste Cayrou * + * copyright (C) 2012-2016 by Jean-Baptiste Cayrou * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -90,14 +90,10 @@ public: virtual std::string classname(void) const { return "GTestModelData"; } virtual std::string type(void) const {return "=== GTestModelData ===";} virtual bool is_constant(void) const {return true;} - virtual double eval(const GEvent& event, - const GObservation& obs) const { - double result = m_modelTps->eval(event.time()); - return result; - } - virtual double eval_gradients(const GEvent& event, - const GObservation& obs) const { - double result = m_modelTps->eval_gradients(event.time()); + virtual double eval(const GEvent& event, + const GObservation& obs, + const bool& gradients ) const { + double result = m_modelTps->eval(event.time(), gradients); return result; } virtual double npred(const GEnergy& obsEng, const GTime& obsTime, -- libgit2 0.21.4