GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelData.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelData.cpp - Abstract virtual data model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-2020 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GModelData.cpp
23  * @brief GModelData class implementation.
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GException.hpp"
33 #include "GModelData.hpp"
34 #include "GVector.hpp"
35 #include "GMatrixSparse.hpp"
36 #include "GObservation.hpp"
37 
38 /* __ Method name definitions ____________________________________________ */
39 #define G_EVAL "GModel::eval(GObservation&, GMatrixSparse*)"
40 
41 /* __ Macros _____________________________________________________________ */
42 
43 /* __ Coding definitions _________________________________________________ */
44 
45 /* __ Debug definitions __________________________________________________ */
46 
47 
48 /*==========================================================================
49  = =
50  = Constructors/destructors =
51  = =
52  ==========================================================================*/
53 
54 /***********************************************************************//**
55  * @brief Void constructor
56  ***************************************************************************/
58 {
59  // Initialise members
60  init_members();
61 
62  // Return
63  return;
64 }
65 
66 
67 /***********************************************************************//**
68  * @brief XML constructor
69  *
70  * @param[in] xml XML element.
71  ***************************************************************************/
73 {
74  // Initialise members
75  init_members();
76 
77  // Return
78  return;
79 }
80 
81 
82 /***********************************************************************//**
83  * @brief Copy constructor
84  *
85  * @param[in] model Data model.
86  ***************************************************************************/
88 {
89  // Initialise members
90  init_members();
91 
92  // Copy members
93  copy_members(model);
94 
95  // Return
96  return;
97 }
98 
99 
100 /***********************************************************************//**
101  * @brief Destructor
102  ***************************************************************************/
104 {
105  // Free members
106  free_members();
107 
108  // Return
109  return;
110 }
111 
112 
113 /*==========================================================================
114  = =
115  = Operators =
116  = =
117  ==========================================================================*/
118 
119 /***********************************************************************//**
120  * @brief Assignment operator
121  *
122  * @param[in] model Data model.
123  * @return Data model.
124  ***************************************************************************/
126 {
127  // Execute only if object is not identical
128  if (this != &model) {
129 
130  // Copy base class members
131  this->GModel::operator=(model);
132 
133  // Free members
134  free_members();
135 
136  // Initialise members
137  init_members();
138 
139  // Copy members
140  copy_members(model);
141 
142  } // endif: object was not identical
143 
144  // Return
145  return *this;
146 }
147 
148 
149 /*==========================================================================
150  = =
151  = Public methods =
152  = =
153  ==========================================================================*/
154 
155 /***********************************************************************//**
156  * @brief Return model values and gradients
157  *
158  * @param[in] obs Observation.
159  * @param[out] gradients Pointer to matrix of gradients.
160  * @return Model values.
161  *
162  * Evaluates the model values and parameter gradients for all events in an
163  * observation. Gradients are only returned if the @p gradients pointer is
164  * not NULL.
165  *
166  * The matrix of gradients is a sparse matrix where the number of rows
167  * corresponds to the number of model parameters and the number of columns
168  * corresponds to the number of events. An exception is thrown if the
169  * dimension of the @p gradients matrix is not compatible with the model
170  * and the observations.
171  ***************************************************************************/
173  GMatrixSparse* gradients) const
174 {
175  // Get number of model parameters and number of events
176  int npars = size();
177  int nevents = obs.events()->size();
178 
179  // Initialise gradients flag
180  bool grad = ((gradients != NULL) && (npars > 0));
181 
182  // Check matrix consistency
183  if (grad) {
184  if (gradients->columns() != npars) {
185  std::string msg = "Number of "+gammalib::str(gradients->columns())+
186  " columns in gradient matrix differs from number "
187  "of "+gammalib::str(npars)+" parameters "
188  "in model. Please provide a compatible gradient "
189  "matrix.";
191  }
192  if (gradients->rows() != nevents) {
193  std::string msg = "Number of "+gammalib::str(gradients->rows())+
194  " rows in gradient matrix differs from number "
195  "of "+gammalib::str(nevents)+" events in "
196  "observation. Please provide a compatible "
197  "gradient matrix.";
199  }
200  }
201 
202  // Initialise values
203  GVector values(nevents);
204 
205  // Initialise temporary vectors to hold gradients
206  GVector* tmp_gradients = NULL;
207  if (grad) {
208 
209  // Allocate temporary vectors to hold the gradients
210  tmp_gradients = new GVector[npars];
211  for (int i = 0; i < npars; ++i) {
212  tmp_gradients[i] = GVector(nevents);
213  }
214 
215  // Signal all parameters that will have analytical gradients. These
216  // are all parameters that are free and for which the model provides
217  // analytical gradients.
218  for (int i = 0; i < npars; ++i) {
219  const GModelPar& par = (*this)[i];
220  if (par.is_free() && par.has_grad()) {
221  obs.computed_gradient(*this, par);
222  }
223  }
224 
225  } // endif: gradients were requested
226 
227  // Loop over events
228  for (int k = 0; k < nevents; ++k) {
229 
230  // Get reference to event
231  const GEvent& event = *((*obs.events())[k]);
232 
233  // Evaluate probability
234  values[k] = eval(event, obs, grad);
235 
236  // Set gradients
237  if (grad) {
238 
239  // Extract relevant parameter gradients
240  for (int i = 0; i < npars; ++i) {
241  const GModelPar& par = (*this)[i];
242  if (par.is_free() && par.has_grad()) {
243  tmp_gradients[i][k] = par.factor_gradient();
244  }
245  }
246 
247  } // endif: looped over gradients
248 
249  } // endfor: looped over events
250 
251  // Post-process gradients
252  if (grad) {
253 
254  // Fill gradients into matrix
255  for (int i = 0; i < npars; ++i) {
256  gradients->column(i, tmp_gradients[i]);
257  }
258 
259  // Delete temporal gradients
260  delete [] tmp_gradients;
261 
262  } // endif: post-processed gradients
263 
264  // Return values
265  return values;
266 }
267 
268 
269 /*==========================================================================
270  = =
271  = Private methods =
272  = =
273  ==========================================================================*/
274 
275 /***********************************************************************//**
276  * @brief Initialise class members
277  ***************************************************************************/
279 {
280  // Return
281  return;
282 }
283 
284 
285 /***********************************************************************//**
286  * @brief Copy class members
287  *
288  * @param[in] model Data model.
289  ***************************************************************************/
291 {
292  // Return
293  return;
294 }
295 
296 
297 /***********************************************************************//**
298  * @brief Delete class members
299  ***************************************************************************/
301 {
302  // Return
303  return;
304 }
Abstract model class.
Definition: GModel.hpp:100
void computed_gradient(const GModel &model, const GModelPar &par) const
Signals that an analytical gradient was computed for a model parameter.
const double & factor_gradient(void) const
Return parameter factor gradient.
Sparse matrix class interface definition.
Abstract interface for the event classes.
Definition: GEvent.hpp:71
XML element node class.
Definition: GXmlElement.hpp:48
void free_members(void)
Delete class members.
Definition: GModelData.cpp:300
#define G_EVAL
Definition: GModelData.cpp:39
Gammalib tools definition.
bool is_free(void) const
Signal if parameter is free.
int size(void) const
Return number of parameters in model.
Definition: GModel.hpp:233
Model parameter class.
Definition: GModelPar.hpp:87
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
const int & rows(void) const
Return number of matrix rows.
Abstract data model base class interface definition.
GModelData(void)
Void constructor.
Definition: GModelData.cpp:57
Abstract data model class.
Definition: GModelData.hpp:55
virtual GModelData & operator=(const GModelData &model)
Assignment operator.
Definition: GModelData.cpp:125
Abstract observation base class.
Vector class interface definition.
Abstract observation base class interface definition.
virtual ~GModelData(void)
Destructor.
Definition: GModelData.cpp:103
void copy_members(const GModelData &model)
Copy class members.
Definition: GModelData.cpp:290
virtual GVector column(const int &column) const
Extract column as vector from matrix.
Sparse matrix class definition.
virtual GEvents * events(void)
Return events.
Exception handler interface definition.
virtual GModel & operator=(const GModel &model)
Assignment operator.
Definition: GModel.cpp:132
Vector class.
Definition: GVector.hpp:46
virtual double eval(const GEvent &event, const GObservation &obs, const bool &gradients=false) const =0
const int & columns(void) const
Return number of matrix columns.
void init_members(void)
Initialise class members.
Definition: GModelData.cpp:278
virtual int size(void) const =0
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489