From 473731bc7db89910e734333853c8905ce12f2f58 Mon Sep 17 00:00:00 2001
From: Michael Mayer <michael.mayer@desy.de>
Date: Tue, 30 Jun 2015 13:12:54 +0200
Subject: [PATCH] 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)

---
 scripts/csspec.py         |  2 +-
 src/cterror/cterror.cpp   |  2 ++
 src/ctulimit/ctulimit.cpp | 59 ++++++++++++++++++++++++++++++++++++++++++-----------------
 src/ctulimit/ctulimit.hpp |  1 +
 test/test_cscripts.sh.in  |  1 +
 5 files changed, 47 insertions(+), 18 deletions(-)

diff --git a/scripts/csspec.py b/scripts/csspec.py
index 5386d74..d5379a9 100755
--- a/scripts/csspec.py
+++ b/scripts/csspec.py
@@ -21,7 +21,7 @@
 import gammalib
 import ctools
 import sys
-import tempfile
+
 
 # ============ #
 # csspec class #
diff --git a/src/cterror/cterror.cpp b/src/cterror/cterror.cpp
index c014762..89583d7 100644
--- a/src/cterror/cterror.cpp
+++ b/src/cterror/cterror.cpp
@@ -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& app)
     m_dlogL     = app.m_dlogL;
     m_best_logL = app.m_best_logL;
     m_model_par = NULL;
+    m_opt = NULL;
 
     // Return
     return;
diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp
index 44c1b6e..92b8110 100644
--- a/src/ctulimit/ctulimit.cpp
+++ b/src/ctulimit/ctulimit.cpp
@@ -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& 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& 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
diff --git a/src/ctulimit/ctulimit.hpp b/src/ctulimit/ctulimit.hpp
index 8eb6727..c1728db 100644
--- a/src/ctulimit/ctulimit.hpp
+++ b/src/ctulimit/ctulimit.hpp
@@ -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
 };
 
 
diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in
index 050aa2f..2010c40 100755
--- a/test/test_cscripts.sh.in
+++ b/test/test_cscripts.sh.in
@@ -44,6 +44,7 @@ rm -rf *.fits *.xml *.dat
 # ========================
 rm -rf pfiles
 mkdir -p pfiles
+cp -r ../src/*/*.par pfiles/
 export PFILES=pfiles
 
 
--
libgit2 0.21.4