From 65cfe455aae323a274ef6e504e7c452076039bbc Mon Sep 17 00:00:00 2001
From: astrojhgu <astrojhgu@ed2142bd-67ad-457f-ba7c-d818d4011675>
Date: Sun, 1 May 2011 17:03:12 +0000
Subject: git-svn-id: file:///home/svn/opt_utilities@191
 ed2142bd-67ad-457f-ba7c-d818d4011675

---
 error_estimator/error_estimator.hpp | 38 ++++++++++++++++++++++++++++++++++++-
 example/test_fitter.cpp             |  3 +--
 interface/type_depository.hpp       |  9 ++++++---
 3 files changed, 44 insertions(+), 6 deletions(-)

diff --git a/error_estimator/error_estimator.hpp b/error_estimator/error_estimator.hpp
index e76fb12..e7a0cc0 100644
--- a/error_estimator/error_estimator.hpp
+++ b/error_estimator/error_estimator.hpp
@@ -8,12 +8,39 @@
 #define ERROR_EST
 #include <core/fitter.hpp>
 #include <core/freeze_param.hpp>
+#include <interface/type_depository.hpp>
+#include <math/num_diff.hpp>
 #include <iostream>
 #include <vector>
 #include <string>
 
 namespace opt_utilities
 {
+  template <typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>
+  typename element_type_trait<Tp>::element_type estimate_error_hessian(fitter<Ty,Tx,Tp,Ts>& fit,const Tstr& pname,const Ts& dchi)
+  {
+    size_t porder=fit.get_param_order(pname);
+    holder<param_modifier<Ty,Tx,Tp,Tstr> > h;
+    try
+      {
+	h.reset(fit.get_param_modifier().clone());
+      }
+    catch(const param_modifier_undefined&)
+      {
+	h.reset(0);
+      }
+    fit.clear_param_modifier();
+    Tp p=fit.get_all_params();
+    Ts a=hessian(fit.get_statistic(),p,porder,porder)/2;
+    if(h.get())
+      {
+	fit.set_param_modifier(*(h.get()));
+      }
+    fit.fit();
+    return std::sqrt(dchi/a);
+  }
+
+
   /**
      \brief calculate the error boundary of a fit, according to the given delta statistic.
      \param fit the fitter that has a sucessful fit result
@@ -24,7 +51,7 @@ namespace opt_utilities
      \param precision determine how precise the error bounds should be determined
    */
   template <typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>
-  void estimate_error(fitter<Ty,Tx,Tp,Ts>& fit,const Tstr& pname,typename element_type_trait<Tp>::element_type& lower,typename element_type_trait<Tp>::element_type& upper,const Ts& dchi,const Ts& precision)
+  void estimate_error_directly(fitter<Ty,Tx,Tp,Ts>& fit,const Tstr& pname,typename element_type_trait<Tp>::element_type& lower,typename element_type_trait<Tp>::element_type& upper,const Ts& dchi,const Ts& precision)
   {
     typedef typename element_type_trait<Tp>::element_type Tpe;
     std::vector<Tstr> pnames;
@@ -175,6 +202,15 @@ namespace opt_utilities
     
     fit.fit();
   }
+
+  template <typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>
+  void estimate_error(fitter<Ty,Tx,Tp,Ts>& fit,const Tstr& pname,typename element_type_trait<Tp>::element_type& lower,typename element_type_trait<Tp>::element_type& upper,const Ts& dchi,const Ts& precision)
+  {
+    typename element_type_trait<Tp>::element_type e=estimate_error_hessian(fit,pname,dchi);
+    lower=fit.get_param_value(pname)-e*2;
+    upper=fit.get_param_value(pname)+e*2;
+    estimate_error_directly(fit,pname,lower,upper,dchi,precision);
+  }
 }
 
 
diff --git a/example/test_fitter.cpp b/example/test_fitter.cpp
index 718d029..47a81ec 100644
--- a/example/test_fitter.cpp
+++ b/example/test_fitter.cpp
@@ -60,7 +60,6 @@ int main()
   
   double lower=3;
   double upper=6;
-  
-  estimate_error(f,std::string("b"),lower,upper,1.,1E-10);
+  estimate_error(f,std::string("k"),lower,upper,1.,1E-10);
   std::cout<<lower<<"\t"<<upper<<endl;
 }
diff --git a/interface/type_depository.hpp b/interface/type_depository.hpp
index d0e0bfb..9bacc9b 100644
--- a/interface/type_depository.hpp
+++ b/interface/type_depository.hpp
@@ -53,7 +53,7 @@ namespace opt_utilities
 
     ~holder()
     {
-      delete ptr;
+      destroy();
     }
 
     holder& operator=(const holder& rhs)
@@ -62,7 +62,7 @@ namespace opt_utilities
 	{
 	  return *this;
 	}
-      delete ptr;
+      destroy();
       ptr=rhs.ptr;
       const_cast<holder&>(rhs).ptr=0;
     }
@@ -77,7 +77,10 @@ namespace opt_utilities
 
     void destroy()
     {
-      delete ptr;
+      if(ptr)
+	{
+	  ptr->destroy();
+	}
     }
 
     T* get()const
-- 
cgit v1.2.2