ctulimit now re-optimises other parameters. I also removed an unused import from…

… csspec and complemented the copy constructor of cterror (the optimiser was not copied and intialised)

ctulimit now re-optimises other parameters. I also removed an unused import from…
… csspec and complemented the copy constructor of cterror (the optimiser was not copied and intialised)
Michael Mayer
1 parent aabec45b
Showing 5 changed files with 47 additions and 18 deletions
scripts/csspec.py
src/cterror/cterror.cpp
src/ctulimit/ctulimit.cpp
src/ctulimit/ctulimit.hpp
test/test_cscripts.sh.in
@@ -21,7 +21,7 @@
 import gammalib
 import ctools
 import sys
-import tempfile
+
  
 # ============ #
 # csspec class #
@@ -412,6 +412,7 @@ void cterror::init_members(void)
     m_dlogL        = 0.0;
     m_best_logL    = 0.0;
     m_model_par    = NULL;
+    m_opt          = NULL;
  
     // Return
     return;
@@ -436,6 +437,7 @@ void cterror::copy_members(const cterror&amp; app)
     m_dlogL     = app.m_dlogL;
     m_best_logL = app.m_best_logL;
     m_model_par = NULL;
+    m_opt = NULL;
  
     // Return
     return;
@@ -232,6 +232,12 @@ void ctulimit::run(void)
     // Save original models
     GModels models_orig = m_obs.models();
  
+    // Allocate a new optimizer
+    m_opt = new GOptimizerLM();
+
+    // Initialise best fit optimizer
+    GOptimizerLM best_opt;
+
     // Save original log-likelihood. If the value is zero it has never been
     // computed hence we compute it now. 
     m_best_logL = m_obs.logL();
@@ -243,24 +249,32 @@ void ctulimit::run(void)
             log.header1("Compute best-fit likelihood");
         }
  
-        // Reoptimize if likelihood was not given before
-        GOptimizerLM* opt = new GOptimizerLM();
-        m_obs.optimize(*opt);
-        m_obs.errors(*opt);
-        m_best_logL = m_obs.logL();
+         // Optimize and save best log-likelihood
+         m_obs.optimize(*m_opt);
+         m_obs.errors(*m_opt);
+         m_best_logL = m_obs.logL();
  
-        // Write optimised model into logger
-        if (logTerse()) {
-            log << m_obs.models() << std::endl;
-        }
+         // Store optimizer for later recovery
+         best_opt = *m_opt;
+
+         // Write optimised model into logger
+         if (logTerse()) {
+             log << *m_opt << std::endl;
+             log << gammalib::parformat("Maximum log likelihood");
+             log << gammalib::str(m_best_logL,3) << std::endl;
+             log << m_obs.models() << std::endl;
+         }
  
     } // endif: likelihood was zero
  
+    // Extract current value
+    double value = m_model_par->factor_value();
+
     // Compute parameter bracketing
-    double value  = m_model_par->value();
-    double error  = m_model_par->error();
-    double parmin = value + m_sigma_min * error;
-    double parmax = value + m_sigma_max * error;
+    double parmin = std::max(m_model_par->factor_min(),
+                             value - m_sigma_min * m_model_par->factor_error());
+    double parmax = std::min(m_model_par->factor_max(),
+                             value + m_sigma_max * m_model_par->factor_error());
  
     // Write header
     if (logTerse()) {
@@ -329,6 +343,9 @@ void ctulimit::run(void)
     // Recover original models
     m_obs.models(models_orig);
  
+    // Recover optimizer
+    *m_opt = best_opt;
+
     // Return
     return;
 }
@@ -401,6 +418,7 @@ void ctulimit::init_members(void)
     m_flux_ulimit  = 0.0;
     m_diff_ulimit  = 0.0;
     m_eflux_ulimit = 0.0;
+    m_opt = NULL;
  
     // Return
     return;
@@ -433,6 +451,7 @@ void ctulimit::copy_members(const ctulimit&amp; app)
     m_diff_ulimit  = app.m_diff_ulimit;
     m_flux_ulimit  = app.m_flux_ulimit;
     m_eflux_ulimit = app.m_eflux_ulimit;
+    m_opt = NULL;
  
     // Extract model parameter
     get_model_parameter();
@@ -662,17 +681,23 @@ double ctulimit::evaluate(const double&amp; value)
     double logL = 0.0;
  
     // Check if given parameter is within boundaries
-    if (value > m_model_par->min() && value < m_model_par->max()) {
+    if (value > m_model_par->factor_min() && value < m_model_par->factor_max()) {
  
         // Change parameter factor
-        m_model_par->value(value);
+        m_model_par->factor_value(value);
  
-        // Evaluate likelihood for new model container
-        m_obs.eval();
+        // Fix parameter
+        m_model_par->fix();
+
+        // Re-optimize
+        m_obs.optimize(*m_opt);
  
         // Retrieve likelihood
         logL = m_obs.logL();
  
+        // Free parameter
+        m_model_par->free();
+
     } // endif: value was inside allowed range
  
     // ... otherwise signal that the parameter went outside the boundaries
@@ -102,6 +102,7 @@ protected:
     double        m_diff_ulimit;  //!< Differential upper limit value
     double        m_flux_ulimit;  //!< Flux upper limit value
     double        m_eflux_ulimit; //!< Energy flux upper limits
+    GOptimizerLM* m_opt;          //!< Optimizer
 };
  
  
@@ -44,6 +44,7 @@ rm -rf *.fits *.xml *.dat
 # ========================
 rm -rf pfiles
 mkdir -p pfiles
+cp -r ../src/*/*.par pfiles/
 export PFILES=pfiles