00001 /*************************************************************************** 00002 * GObservations.cpp - Observation container class * 00003 * ----------------------------------------------------------------------- * 00004 * copyright (C) 2009-2016 by Juergen Knoedlseder * 00005 * ----------------------------------------------------------------------- * 00006 * * 00007 * This program is free software: you can redistribute it and/or modify * 00008 * it under the terms of the GNU General Public License as published by * 00009 * the Free Software Foundation, either version 3 of the License, or * 00010 * (at your option) any later version. * 00011 * * 00012 * This program is distributed in the hope that it will be useful, * 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 00015 * GNU General Public License for more details. * 00016 * * 00017 * You should have received a copy of the GNU General Public License * 00018 * along with this program. If not, see <http://www.gnu.org/licenses/>. * 00019 * * 00020 ***************************************************************************/ 00021 /** 00022 * @file GObservations.cpp 00023 * @brief Observations container class implementation 00024 * @author Juergen Knoedlseder 00025 */ 00026 00027 /* __ Includes ___________________________________________________________ */ 00028 #ifdef HAVE_CONFIG_H 00029 #include <config.h> 00030 #endif 00031 #include "GTools.hpp" 00032 #include "GException.hpp" 00033 #include "GFilename.hpp" 00034 #include "GObservations.hpp" 00035 #include "GObservationRegistry.hpp" 00036 #include "GMatrixSparse.hpp" 00037 #include "GOptimizer.hpp" 00038 00039 /* __ Method name definitions ____________________________________________ */ 00040 #define G_AT "GObservations::at(int&)" 00041 #define G_SET "GObservations::set(int&, GObservation&)" 00042 #define G_APPEND "GObservations::append(GObservation&)" 00043 #define G_INSERT "GObservations::insert(int&, GObservation&)" 00044 #define G_REMOVE "GObservations::remove(int&)" 00045 #define G_EXTEND "GObservations::extend(GObservations&)" 00046 #define G_READ "GObservations::read(GXml&)" 00047 00048 /* __ Macros _____________________________________________________________ */ 00049 00050 /* __ Coding definitions _________________________________________________ */ 00051 00052 /* __ Debug definitions __________________________________________________ */ 00053 00054 00055 /*========================================================================== 00056 = = 00057 = Constructors/destructors = 00058 = = 00059 ==========================================================================*/ 00060 00061 /***********************************************************************//** 00062 * @brief Void constructor 00063 ***************************************************************************/ 00064 GObservations::GObservations(void) 00065 { 00066 // Initialise members 00067 init_members(); 00068 00069 // Return 00070 return; 00071 } 00072 00073 00074 /***********************************************************************//** 00075 * @brief Copy constructor 00076 * 00077 * @param obs Observation container. 00078 ***************************************************************************/ 00079 GObservations::GObservations(const GObservations& obs) 00080 { 00081 // Initialise members 00082 init_members(); 00083 00084 // Copy members 00085 copy_members(obs); 00086 00087 // Return 00088 return; 00089 } 00090 00091 00092 /***********************************************************************//** 00093 * @brief Load constructor 00094 * 00095 * @param[in] filename XML filename. 00096 * 00097 * Construct the observation container by loading all relevant information 00098 * from a XML file. Please refer to the read() method for more information 00099 * about the structure of the XML file. 00100 ***************************************************************************/ 00101 GObservations::GObservations(const GFilename& filename) 00102 { 00103 // Initialise members 00104 init_members(); 00105 00106 // Load XML file 00107 load(filename); 00108 00109 // Return 00110 return; 00111 } 00112 00113 00114 /***********************************************************************//** 00115 * @brief Destructor 00116 ***************************************************************************/ 00117 GObservations::~GObservations(void) 00118 { 00119 // Free members 00120 free_members(); 00121 00122 // Return 00123 return; 00124 } 00125 00126 00127 /*========================================================================== 00128 = = 00129 = Operators = 00130 = = 00131 ==========================================================================*/ 00132 00133 /***********************************************************************//** 00134 * @brief Assignment operator 00135 * 00136 * @param[in] obs Observation container. 00137 * @return Observation container. 00138 ***************************************************************************/ 00139 GObservations& GObservations::operator=(const GObservations& obs) 00140 { 00141 // Execute only if object is not identical 00142 if (this != &obs) { 00143 00144 // Free members 00145 free_members(); 00146 00147 // Initialise members 00148 init_members(); 00149 00150 // Copy members 00151 copy_members(obs); 00152 00153 } // endif: object was not identical 00154 00155 // Return this object 00156 return *this; 00157 } 00158 00159 00160 /*========================================================================== 00161 = = 00162 = Public methods = 00163 = = 00164 ==========================================================================*/ 00165 00166 /***********************************************************************//** 00167 * @brief Clear observations 00168 ***************************************************************************/ 00169 void GObservations::clear(void) 00170 { 00171 // Free members 00172 free_members(); 00173 00174 // Initialise members 00175 init_members(); 00176 00177 // Return 00178 return; 00179 } 00180 00181 00182 /***********************************************************************//** 00183 * @brief Clone observations 00184 * 00185 * @return Pointer to deep copy of observation container. 00186 ***************************************************************************/ 00187 GObservations* GObservations::clone(void) const 00188 { 00189 // Clone observations 00190 return new GObservations(*this); 00191 } 00192 00193 00194 /***********************************************************************//** 00195 * @brief Return pointer to observation 00196 * 00197 * @param[in] index Observation index [0,...,size()-1]. 00198 * @return Observation. 00199 * 00200 * @exception GException::out_of_range 00201 * Operation index is out of range. 00202 * 00203 * Returns a pointer to the observation with specified @p index. 00204 ***************************************************************************/ 00205 GObservation* GObservations::at(const int& index) 00206 { 00207 // Raise exception if index is out of range 00208 if (index < 0 || index >= size()) { 00209 throw GException::out_of_range(G_AT, index, 0, size()-1); 00210 } 00211 00212 // Return pointer 00213 return m_obs[index]; 00214 } 00215 00216 00217 /***********************************************************************//** 00218 * @brief Return pointer to observation (const version) 00219 * 00220 * @param[in] index Observation index [0,...,size()-1]. 00221 * 00222 * @exception GException::out_of_range 00223 * Operation index is out of range. 00224 * 00225 * Returns a const pointer to the observation with specified @p index. 00226 ***************************************************************************/ 00227 const GObservation* GObservations::at(const int& index) const 00228 { 00229 // Raise exception if index is out of range 00230 if (index < 0 || index >= size()) { 00231 throw GException::out_of_range(G_AT, index, 0, size()-1); 00232 } 00233 00234 // Return pointer 00235 return m_obs[index]; 00236 } 00237 00238 00239 /***********************************************************************//** 00240 * @brief Set observation in container 00241 * 00242 * @param[in] index Observation index [0,...,size()-1]. 00243 * @param[in] obs Observation. 00244 * @return Pointer to deep copy of observation. 00245 * 00246 * @exception GException::out_of_range 00247 * Observation index is out of range. 00248 * @exception GException::invalid_value 00249 * Observation with same instrument and identifier already 00250 * exists in container. 00251 * 00252 * Set a deep copy and observation @p obs at the specified @p index in the 00253 * container. 00254 ***************************************************************************/ 00255 GObservation* GObservations::set(const int& index, const GObservation& obs) 00256 { 00257 // Compile option: raise exception if index is out of range 00258 #if defined(G_RANGE_CHECK) 00259 if (index < 0 || index >= size()) { 00260 throw GException::out_of_range(G_SET, index, 0, size()-1); 00261 } 00262 #endif 00263 00264 // Raise an exception if an observation with specified instrument and 00265 // identifier already exists 00266 int inx = get_index(obs.instrument(), obs.id()); 00267 if (inx != -1 && inx != index) { 00268 std::string msg = 00269 "Attempt to set \""+obs.instrument()+"\" observation with" 00270 " identifier \""+obs.id()+"\" in observation container at" 00271 " index "+gammalib::str(index)+", but an observation with the" 00272 " same attributes exists already at index "+gammalib::str(inx)+ 00273 " in the container.\n" 00274 "Every observation for a given instrument in the observation" 00275 " container needs a unique identifier."; 00276 throw GException::invalid_value(G_SET, msg); 00277 } 00278 00279 // Delete existing observation 00280 if (m_obs[index] != NULL) delete m_obs[index]; 00281 00282 // Store pointer to a deep copy of the observation 00283 m_obs[index] = obs.clone(); 00284 00285 // Return pointer 00286 return m_obs[index]; 00287 } 00288 00289 00290 /***********************************************************************//** 00291 * @brief Append observation to container 00292 * 00293 * @param[in] obs Observation. 00294 * @return Pointer to deep copy of observation. 00295 * 00296 * @exception GException::invalid_value 00297 * Observation with same instrument and identifier already 00298 * exists in container. 00299 * 00300 * Appends a deep copy of an observation to the container. 00301 ***************************************************************************/ 00302 GObservation* GObservations::append(const GObservation& obs) 00303 { 00304 // Raise an exception if an observation with specified instrument and 00305 // identifier already exists 00306 int inx = get_index(obs.instrument(), obs.id()); 00307 if (inx != -1) { 00308 std::string msg = 00309 "Attempt to append \""+obs.instrument()+"\" observation with" 00310 " identifier \""+obs.id()+"\" to observation container, but an" 00311 " observation with the same attributes exists already at" 00312 " index "+gammalib::str(inx)+" in the container.\n" 00313 "Every observation for a given instrument in the observation" 00314 " container needs a unique identifier."; 00315 throw GException::invalid_value(G_APPEND, msg); 00316 } 00317 00318 // Clone observation 00319 GObservation* ptr = obs.clone(); 00320 00321 // Append to list 00322 m_obs.push_back(ptr); 00323 00324 // Return pointer 00325 return ptr; 00326 } 00327 00328 00329 /***********************************************************************//** 00330 * @brief Insert observation into container 00331 * 00332 * @param[in] index Observation index [0,...,size()-1]. 00333 * @param[in] obs Observation. 00334 * @return Pointer to deep copy of observation. 00335 * 00336 * @exception GException::out_of_range 00337 * Observation index is out of range. 00338 * @exception GException::invalid_value 00339 * Observation with same instrument and identifier already 00340 * exists in container. 00341 * 00342 * Inserts a deep copy of an observation into the container before the 00343 * observation with the specified @p index. 00344 ***************************************************************************/ 00345 GObservation* GObservations::insert(const int& index, const GObservation& obs) 00346 { 00347 // Compile option: raise exception if index is out of range 00348 #if defined(G_RANGE_CHECK) 00349 if (is_empty()) { 00350 if (index > 0) { 00351 throw GException::out_of_range(G_INSERT, index, 0, size()-1); 00352 } 00353 } 00354 else { 00355 if (index < 0 || index >= size()) { 00356 throw GException::out_of_range(G_INSERT, index, 0, size()-1); 00357 } 00358 } 00359 #endif 00360 00361 // Raise an exception if an observation with specified instrument and 00362 // identifier already exists 00363 int inx = get_index(obs.instrument(), obs.id()); 00364 if (inx != -1 && inx != index) { 00365 std::string msg = 00366 "Attempt to insert \""+obs.instrument()+"\" observation with" 00367 " identifier \""+obs.id()+"\" in observation container before" 00368 " index "+gammalib::str(index)+", but an observation with the" 00369 " same attributes exists already at index "+gammalib::str(inx)+ 00370 " in the container.\n" 00371 "Every observation for a given instrument in the observation" 00372 " container needs a unique identifier."; 00373 throw GException::invalid_value(G_INSERT, msg); 00374 } 00375 00376 // Clone observation 00377 GObservation* ptr = obs.clone(); 00378 00379 // Clone observation and insert into list 00380 m_obs.insert(m_obs.begin()+index, ptr); 00381 00382 // Return pointer 00383 return ptr; 00384 } 00385 00386 00387 /***********************************************************************//** 00388 * @brief Remove observation from container 00389 * 00390 * @param[in] index Observation index [0,...,size()-1]. 00391 * 00392 * @exception GException::out_of_range 00393 * Observation index is out of range. 00394 * 00395 * Removes observation of specified @p index from the container. 00396 ***************************************************************************/ 00397 void GObservations::remove(const int& index) 00398 { 00399 // Compile option: If index is outside boundary then raise exception 00400 #if defined(G_RANGE_CHECK) 00401 if (index < 0 || index >= size()) { 00402 throw GException::out_of_range(G_REMOVE, index, 0, size()-1); 00403 } 00404 #endif 00405 00406 // Delete observation 00407 delete m_obs[index]; 00408 00409 // Erase observation component from container 00410 m_obs.erase(m_obs.begin() + index); 00411 00412 // Return 00413 return; 00414 } 00415 00416 00417 /***********************************************************************//** 00418 * @brief Append observations from observation container 00419 * 00420 * @param[in] obs Observations. 00421 * 00422 * @exception GException::invalid_value 00423 * Observation with same instrument and identifier already 00424 * exists in container. 00425 * 00426 * Appends deep copies of observations to the container. 00427 ***************************************************************************/ 00428 void GObservations::extend(const GObservations& obs) 00429 { 00430 // Do nothing if observation container is empty 00431 if (!obs.is_empty()) { 00432 00433 // Get size. Note that we extract the size first to avoid an 00434 // endless loop that arises when a container is appended to 00435 // itself. 00436 int num = obs.size(); 00437 00438 // Reserve enough space 00439 reserve(size() + num); 00440 00441 // Loop over all observations, clone them and append them to the 00442 // list 00443 for (int i = 0; i < num; ++i) { 00444 00445 // Raise an exception if an observation with specified 00446 // instrument and identifier already exists 00447 int inx = get_index(obs[i]->instrument(), obs[i]->id()); 00448 if (inx != -1) { 00449 std::string msg = 00450 "Attempt to append \""+obs[i]->instrument()+"\"" 00451 " observation with identifier \""+obs[i]->id()+"\"" 00452 " to observation container, but an observation with the" 00453 " same attributes exists already at index "+ 00454 gammalib::str(inx)+" in the container.\n" 00455 "Every observation for a given instrument in the" 00456 " observation container needs a unique identifier."; 00457 throw GException::invalid_value(G_EXTEND, msg); 00458 } 00459 00460 // Append observation to container 00461 m_obs.push_back(obs[i]->clone()); 00462 } 00463 00464 } // endif: observation container was not empty 00465 00466 // Return 00467 return; 00468 } 00469 00470 00471 /***********************************************************************//** 00472 * @brief Signals if observation exists 00473 * 00474 * @param[in] instrument Instrument. 00475 * @param[in] id Observation identifier. 00476 * @return True if observation with specified @p instrument and identifier 00477 * @p id exists. 00478 * 00479 * Searches all observations for a match with the specified @p instrument 00480 * and @p identifier. If the specified attributes have been found, true is 00481 * returned. 00482 ***************************************************************************/ 00483 bool GObservations::contains(const std::string& instrument, 00484 const std::string& id) const 00485 { 00486 // Get observation index 00487 int index = get_index(instrument, id); 00488 00489 // Return 00490 return (index != -1); 00491 } 00492 00493 00494 /***********************************************************************//** 00495 * @brief Load observations from XML file 00496 * 00497 * @param[in] filename XML filename. 00498 * 00499 * Loads observation from a XML file into the container. Please refer to the 00500 * read() method for more information about the structure of the XML file. 00501 ***************************************************************************/ 00502 void GObservations::load(const GFilename& filename) 00503 { 00504 // Clear any existing observations 00505 clear(); 00506 00507 // Load XML document 00508 GXml xml(filename.url()); 00509 00510 // Read observations from XML document 00511 read(xml); 00512 00513 // Return 00514 return; 00515 } 00516 00517 00518 /***********************************************************************//** 00519 * @brief Save observations into XML file 00520 * 00521 * @param[in] filename XML filename. 00522 * 00523 * Saves observations into a XML file. Please refer to the read() method for 00524 * more information about the structure of the XML file. 00525 ***************************************************************************/ 00526 void GObservations::save(const GFilename& filename) const 00527 { 00528 // Declare empty XML document 00529 GXml xml; 00530 00531 // Write observations into XML file 00532 write(xml); 00533 00534 // Save XML document 00535 xml.save(filename); 00536 00537 // Return 00538 return; 00539 } 00540 00541 00542 /***********************************************************************//** 00543 * @brief Read observations from XML document 00544 * 00545 * @param[in] xml XML document. 00546 * 00547 * @exception GException::invalid_instrument 00548 * Invalid instrument encountered in XML file. 00549 * 00550 * Reads observations from the first observation list that is found in the 00551 * XML document. The decoding of the instrument specific observation 00552 * definition is done within the observation's GObservation::read() method. 00553 * The following file structure is expected: 00554 * 00555 * <observation_list title="observation library"> 00556 * <observation name="..." id="..." instrument="..."> 00557 * ... 00558 * </observation> 00559 * <observation name="..." id="..." instrument="..."> 00560 * ... 00561 * </observation> 00562 * ... 00563 * </observation_list> 00564 * 00565 * The @p name and @p id attributes allow for a unique identification of an 00566 * observation within the observation container. The @p instrument 00567 * attributes specifies the instrument to which the observation applies. 00568 * This attribute will be used to allocate the appropriate instrument 00569 * specific GObservation class variant by using the GObservationRegistry 00570 * class. 00571 * 00572 * The structure within the @p observation tag is defined by the instrument 00573 * specific GObservation class. 00574 * 00575 * @todo Observation names and IDs are not verified so far for uniqueness. 00576 * This would be required to achieve an unambiguous update of parameters 00577 * in an already existing XML file when using the write method. 00578 ***************************************************************************/ 00579 void GObservations::read(const GXml& xml) 00580 { 00581 // Get pointer on observation library 00582 const GXmlElement* lib = xml.element("observation_list", 0); 00583 00584 // Loop over all observations 00585 int n = lib->elements("observation"); 00586 for (int i = 0; i < n; ++i) { 00587 00588 // Get pointer on observation 00589 const GXmlElement* obs = lib->element("observation", i); 00590 00591 // Get attributes 00592 std::string name = obs->attribute("name"); 00593 std::string id = obs->attribute("id"); 00594 std::string instrument = obs->attribute("instrument"); 00595 00596 // Allocate instrument specific observation 00597 GObservationRegistry registry; 00598 GObservation* ptr = registry.alloc(instrument); 00599 00600 // If observation is valid then read its definition from XML file 00601 if (ptr != NULL) { 00602 00603 // Read definition 00604 ptr->read(*obs); 00605 00606 // Set attributes 00607 ptr->name(name); 00608 ptr->id(id); 00609 00610 } // endif: observation was valid 00611 00612 // ... otherwise throw an exception 00613 else { 00614 throw GException::invalid_instrument(G_READ, instrument); 00615 } 00616 00617 // Append observation to container 00618 append(*ptr); 00619 00620 // Free observation (the append method made a deep copy) 00621 delete ptr; 00622 00623 } // endfor: looped over all observations 00624 00625 // Return 00626 return; 00627 } 00628 00629 00630 /***********************************************************************//** 00631 * @brief Write observations into XML document 00632 * 00633 * @param[in] xml XML document. 00634 * 00635 * Write observations into the first observation library that is found in the 00636 * XML document. In case that no observation library exists, one is added to 00637 * the document. Please refer to the read() method for more information 00638 * about the structure of the XML document. 00639 ***************************************************************************/ 00640 void GObservations::write(GXml& xml) const 00641 { 00642 // If there is no observation library then append one 00643 if (xml.elements("observation_list") == 0) { 00644 xml.append(GXmlElement("observation_list title=\"observation list\"")); 00645 } 00646 00647 // Get pointer on observation library 00648 GXmlElement* lib = xml.element("observation_list", 0); 00649 00650 // Write all observations into library 00651 for (int i = 0; i < size(); ++i) { 00652 00653 // Initialise pointer on observation 00654 GXmlElement* obs = NULL; 00655 00656 // Search corresponding observation 00657 int n = xml.elements("observation"); 00658 for (int k = 0; k < n; ++k) { 00659 GXmlElement* element = xml.element("observation", k); 00660 if (element->attribute("name") == m_obs[i]->name() && 00661 element->attribute("id") == m_obs[i]->id() && 00662 element->attribute("instrument") == m_obs[i]->instrument()) { 00663 obs = element; 00664 break; 00665 } 00666 } 00667 00668 // If no observation with corresponding name, ID and instrument was 00669 // found then append one now 00670 if (obs == NULL) { 00671 obs = lib->append("observation"); 00672 obs->attribute("name", m_obs[i]->name()); 00673 obs->attribute("id", m_obs[i]->id()); 00674 obs->attribute("instrument", m_obs[i]->instrument()); 00675 } 00676 00677 // Write now observation 00678 m_obs[i]->write(*obs); 00679 00680 } // endfor: looped over all observaitons 00681 00682 // Return 00683 return; 00684 } 00685 00686 00687 /***********************************************************************//** 00688 * @brief Load models from XML file 00689 * 00690 * @param[in] filename XML filename. 00691 * 00692 * Loads the models from a XML file. Please refer to the GModels::read() 00693 * method for more information about the expected structure of the XML 00694 * file. 00695 ***************************************************************************/ 00696 void GObservations::models(const GFilename& filename) 00697 { 00698 // Load models 00699 m_models.load(filename); 00700 00701 // Return 00702 return; 00703 } 00704 00705 00706 /***********************************************************************//** 00707 * @brief Optimize model parameters using optimizer 00708 * 00709 * @param[in] opt Optimizer. 00710 * 00711 * Optimizes the free parameters of the models by using the optimizer 00712 * that has been provided by the @p opt argument. 00713 ***************************************************************************/ 00714 void GObservations::optimize(GOptimizer& opt) 00715 { 00716 // Extract optimizer parameter container from model container 00717 GOptimizerPars pars = m_models.pars(); 00718 00719 // Optimize model parameters 00720 opt.optimize(m_fct, pars); 00721 00722 // Return 00723 return; 00724 } 00725 00726 00727 /***********************************************************************//** 00728 * @brief Computes parameter errors using optimizer 00729 * 00730 * @param[in] opt Optimizer. 00731 * 00732 * Calculates the errors of the free parameters of the models by using 00733 * the optimizer that has been provided by the @p opt argument. 00734 ***************************************************************************/ 00735 void GObservations::errors(GOptimizer& opt) 00736 { 00737 // Extract optimizer parameter container from model container 00738 GOptimizerPars pars = m_models.pars(); 00739 00740 // Compute model parameter errors 00741 opt.errors(m_fct, pars); 00742 00743 // Return 00744 return; 00745 } 00746 00747 /***********************************************************************//** 00748 * @brief Computes parameter errors using hessian matrix and optimizer 00749 * 00750 * Calculates the errors of the free model parameters as the square roots 00751 * of the diagonal elements of the inverse Hessian matrix. The Hessian 00752 * matrix is computed using finite differences. 00753 ***************************************************************************/ 00754 void GObservations::errors_hessian(void) 00755 { 00756 // Extract optimizer parameter container from model container 00757 GOptimizerPars pars = m_models.pars(); 00758 00759 // Get number of parameters 00760 int npars = pars.size(); 00761 00762 // Compute Hessian matrix 00763 GMatrixSparse hessian = m_fct.hessian(pars); 00764 00765 // Signal no diagonal element loading 00766 bool diag_loaded = false; 00767 00768 // Loop over error computation (maximum 2 passes) 00769 for (int i = 0; i < 2; ++i) { 00770 00771 // Solve: hessian * X = unit 00772 try { 00773 GMatrixSparse decomposition = hessian.cholesky_decompose(true); 00774 GVector unit(npars); 00775 for (int ipar = 0; ipar < npars; ++ipar) { 00776 unit[ipar] = 1.0; 00777 GVector x = decomposition.cholesky_solver(unit, true); 00778 if (x[ipar] >= 0.0) { 00779 pars[ipar]->factor_error(sqrt(x[ipar])); 00780 } 00781 else { 00782 pars[ipar]->factor_error(0.0); 00783 } 00784 unit[ipar] = 0.0; 00785 } 00786 } 00787 catch (GException::matrix_zero &e) { 00788 std::cout << "GObservations::errors_hessian: " 00789 << "All hessian matrix elements are zero." 00790 << std::endl; 00791 break; 00792 } 00793 catch (GException::matrix_not_pos_definite &e) { 00794 00795 // Load diagonal if this has not yet been tried 00796 if (!diag_loaded) { 00797 00798 // Flag errors as inaccurate 00799 std::cout << "Non-Positive definite hessian matrix encountered." 00800 << std::endl; 00801 std::cout << "Load diagonal elements with 1e-10." 00802 << " Fit errors may be inaccurate." 00803 << std::endl; 00804 00805 // Try now with diagonal loaded matrix 00806 for (int ipar = 0; ipar < npars; ++ipar) { 00807 hessian(ipar,ipar) += 1.0e-10; 00808 } 00809 00810 // Signal loading 00811 diag_loaded = true; 00812 00813 // Try again 00814 continue; 00815 00816 } // endif: diagonal has not yet been loaded 00817 00818 // ... otherwise signal an error 00819 else { 00820 std::cout << "Non-Positive definite hessian matrix encountered," 00821 << " even after diagonal loading." << std::endl; 00822 break; 00823 } 00824 } 00825 catch (std::exception &e) { 00826 throw; 00827 } 00828 00829 // If no error occured then break now 00830 break; 00831 00832 } // endfor: looped over error computation 00833 00834 // Return 00835 return; 00836 } 00837 00838 00839 /***********************************************************************//** 00840 * @brief Evaluate function 00841 * 00842 * Evaluates the likelihood funtion at the actual set of parameters. 00843 ***************************************************************************/ 00844 void GObservations::eval(void) 00845 { 00846 // Extract optimizer parameter container from model container 00847 GOptimizerPars pars = m_models.pars(); 00848 00849 // Compute log-likelihood 00850 m_fct.eval(pars); 00851 00852 // Return 00853 return; 00854 } 00855 00856 00857 /***********************************************************************//** 00858 * @brief Return total number of observed events 00859 * 00860 * @return Total number of observed events. 00861 * 00862 * Returns the total number of observed events that is container in the 00863 * observation container. 00864 ***************************************************************************/ 00865 int GObservations::nobserved(void) const 00866 { 00867 // Initialise number of observed events 00868 int nobserved = 0; 00869 00870 // Compute number of observed events 00871 for (int i = 0; i < size(); ++i) { 00872 nobserved += (*this)[i]->nobserved(); 00873 } 00874 00875 // Return number of observed events 00876 return nobserved; 00877 } 00878 00879 00880 /***********************************************************************//** 00881 * @brief Print observation list information 00882 * 00883 * @param[in] chatter Chattiness (defaults to NORMAL). 00884 * @return String containing observation list information 00885 ***************************************************************************/ 00886 std::string GObservations::print(const GChatter& chatter) const 00887 { 00888 // Initialise result string 00889 std::string result; 00890 00891 // Continue only if chatter is not silent 00892 if (chatter != SILENT) { 00893 00894 // Append header 00895 result.append("=== GObservations ==="); 00896 00897 // Append information 00898 result.append("\n"+gammalib::parformat("Number of observations")); 00899 result.append(gammalib::str(size())); 00900 result.append("\n"+gammalib::parformat("Number of models")); 00901 result.append(gammalib::str(m_models.size())); 00902 result.append("\n"+gammalib::parformat("Number of observed events")); 00903 result.append(gammalib::str(nobserved())); 00904 result.append("\n"+gammalib::parformat("Number of predicted events")); 00905 result.append(gammalib::str(npred())); 00906 00907 // EXPLICIT: Append observations 00908 if (chatter >= EXPLICIT) { 00909 for (int i = 0; i < size(); ++i) { 00910 result.append("\n"); 00911 result.append((*this)[i]->print()); 00912 } 00913 } // endif: chatter was explicit 00914 00915 // VERBOSE: Append models 00916 if (chatter == VERBOSE) { 00917 result.append("\n"+m_models.print()); 00918 } // endif: chatter was verbose 00919 00920 } // endif: chatter was not silent 00921 00922 // Return result 00923 return result; 00924 } 00925 00926 00927 /*========================================================================== 00928 = = 00929 = Private methods = 00930 = = 00931 ==========================================================================*/ 00932 00933 /***********************************************************************//** 00934 * @brief Initialise class members 00935 ***************************************************************************/ 00936 void GObservations::init_members(void) 00937 { 00938 // Initialise members 00939 m_obs.clear(); 00940 m_models.clear(); 00941 m_fct.set(this); //!< Makes sure that optimizer points to this instance 00942 00943 // Return 00944 return; 00945 } 00946 00947 00948 /***********************************************************************//** 00949 * @brief Copy class members 00950 * 00951 * @param[in] obs Observation container. 00952 * 00953 * Copies all class members. Deep copies are created for the observations, 00954 * which is what is needed to have a real physical copy of the members. 00955 ***************************************************************************/ 00956 void GObservations::copy_members(const GObservations& obs) 00957 { 00958 // Copy attributes 00959 m_models = obs.m_models; 00960 m_fct = obs.m_fct; 00961 00962 // Copy observations 00963 m_obs.clear(); 00964 for (int i = 0; i < obs.m_obs.size(); ++i) { 00965 m_obs.push_back((obs.m_obs[i]->clone())); 00966 } 00967 00968 // Set likelihood function pointer to this instance. This makes sure 00969 // that the optimizer points to this instance 00970 m_fct.set(this); 00971 00972 // Return 00973 return; 00974 } 00975 00976 00977 /***********************************************************************//** 00978 * @brief Delete class members 00979 * 00980 * Deallocates all observations. Since container classes that hold pointers 00981 * need to handle the proper deallocation of memory, we loop here over all 00982 * pointers and make sure that we deallocate all observations in the 00983 * container. 00984 ***************************************************************************/ 00985 void GObservations::free_members(void) 00986 { 00987 // Free observations 00988 for (int i = 0; i < m_obs.size(); ++i) { 00989 delete m_obs[i]; 00990 m_obs[i] = NULL; 00991 } 00992 00993 // Return 00994 return; 00995 } 00996 00997 00998 /***********************************************************************//** 00999 * @brief Return observation index by instrument and identifier 01000 * 01001 * @param[in] instrument Instrument. 01002 * @param[in] id Observation identifier. 01003 * @return Observation index (-1 of instrument and observation identifier 01004 * has not been found) 01005 * 01006 * Returns observation index based on the specified @p instrument and 01007 * observation identifier @p id. If no observation with the specified 01008 * attributes has been found, the method returns -1. 01009 ***************************************************************************/ 01010 int GObservations::get_index(const std::string& instrument, 01011 const std::string& id) const 01012 { 01013 // Initialise index 01014 int index = -1; 01015 01016 // Search observation with specified instrument and id 01017 for (int i = 0; i < size(); ++i) { 01018 if ((m_obs[i]->instrument() == instrument) && 01019 (m_obs[i]->id() == id)) { 01020 index = i; 01021 break; 01022 } 01023 } 01024 01025 // Return index 01026 return index; 01027 }