From c38a8d7c3cca296c9f9a70bc3504075b416c68b6 Mon Sep 17 00:00:00 2001
From: astrojhgu <astrojhgu@ed2142bd-67ad-457f-ba7c-d818d4011675>
Date: Tue, 4 Jan 2011 14:23:16 +0000
Subject: git-svn-id: file:///home/svn/opt_utilities@152
 ed2142bd-67ad-457f-ba7c-d818d4011675

---
 core/fitter.hpp                     |   9 +++
 core/freeze_param.hpp               |   4 +-
 error_estimator/error_estimator.hpp | 136 ++++++++++++++++++++++++++++++++++++
 3 files changed, 147 insertions(+), 2 deletions(-)
 create mode 100644 error_estimator/error_estimator.hpp

diff --git a/core/fitter.hpp b/core/fitter.hpp
index 9caae67..5ec85c5 100644
--- a/core/fitter.hpp
+++ b/core/fitter.hpp
@@ -1266,6 +1266,15 @@ namespace opt_utilities
 	}
       return *(this->p_statistic);
     }
+
+    /**
+       \return current statstic value
+     */
+    Ts get_statistic_value()
+    {
+      Tp current_params(get_model().get_all_params());
+      return get_statistic().eval(get_model().deform_param(current_params));
+    }
     
     /**
        Get the optimization method that used
diff --git a/core/freeze_param.hpp b/core/freeze_param.hpp
index 015398e..abd8b17 100644
--- a/core/freeze_param.hpp
+++ b/core/freeze_param.hpp
@@ -150,9 +150,9 @@ namespace opt_utilities
     {
       if(param_names.find(name)==param_names.end())
 	{
-	  return "thawed";
+	  return Tstr("thawed");
 	}
-      return "frozen";
+      return Tstr("frozen");
     }
 
   public:
diff --git a/error_estimator/error_estimator.hpp b/error_estimator/error_estimator.hpp
new file mode 100644
index 0000000..1e9fa84
--- /dev/null
+++ b/error_estimator/error_estimator.hpp
@@ -0,0 +1,136 @@
+#ifndef ERROR_EST
+#define ERROR_EST
+#include <core/fitter.hpp>
+#include <iostream>
+
+namespace opt_utilities
+{
+  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)
+  {
+    
+    typedef typename element_type_trait<Tp>::element_type Tpe;
+    fit.fit();
+    if(fit.report_param_status(pname)=="frozen")
+      {
+	return;
+      }
+    try
+      {
+	freeze_param<Ty,Tx,Tp,Tstr>* pfp=dynamic_cast<freeze_param<Ty,Tx,Tp,Tstr>*>(&fit.get_param_modifier());
+	if(pfp==0)
+	  {
+	    assert(0);
+	  }
+	*pfp+=freeze_param<Ty,Tx,Tp,Tstr>(pname);
+      }
+    catch(const param_modifier_undefined&)
+      {
+	fit.set_param_modifier(freeze_param<Ty,Tx,Tp,Tstr>(pname));
+      }
+    
+    Tpe current_value=
+      fit.get_param_value(pname);
+    Tp current_params(fit.get_model().get_all_params());
+    Ts current_chi=fit.get_statistic_value();
+    Tpe upper1=current_value;
+    Tpe step=upper-current_value;
+    while(1)
+      {
+	fit.set_param_value(pname,upper);
+	fit.fit();
+	Ts chi=fit.get_statistic_value();
+	cerr<<upper<<"\t"<<chi<<std::endl;
+	if(chi>current_chi+dchi)
+	  {
+	    break;
+	  }
+	upper+=step;
+      }
+    
+    Ts target_chi=current_chi+dchi;
+    while(upper-upper1>1E-10)
+      {
+	std::cerr<<upper<<"\t"<<upper1<<std::endl;
+	fit.set_param_value(pname,upper1);
+	fit.fit();
+	Ts s1=fit.get_statistic_value();
+	fit.set_param_value(pname,upper);
+	fit.fit();
+	Ts s=fit.get_statistic_value();
+	fit.set_param_value(pname,(upper+upper1)/2);
+	fit.fit();
+	Ts s01=fit.get_statistic_value();
+	if((s1-target_chi)*(s01-target_chi)<0&&(s-target_chi)*(s01-target_chi)>=0)
+	  {
+	    upper=(upper+upper1)/2;
+	  }
+	else if((s1-target_chi)*(s01-target_chi)>0&&(s-target_chi)*(s01-target_chi)<=0)
+	  {
+	    upper1=(upper+upper1)/2;
+	  }
+	else
+	  {
+	    std::cerr<<"strange statistic structure"<<std::endl;
+	    break;
+	  }
+      }
+    
+    fit.set_param_value(pname,current_value);
+    fit.fit();
+    
+    ////
+    Tpe lower1=current_value;
+    step=current_value-lower;
+    while(1)
+      {
+	fit.set_param_value(pname,lower);
+	fit.fit();
+	Ts chi=fit.get_statistic_value();
+	std::cerr<<lower<<"\t"<<chi<<std::endl;
+	if(chi>current_chi+dchi)
+	  {
+	    break;
+	  }
+	lower-=step;
+      }
+    
+    fit.set_param_value(pname,current_value);
+    fit.fit();
+    
+    while(lower1-lower>1E-10)
+      {
+	fit.set_param_value(pname,lower1);
+	fit.fit();
+	Ts s1=fit.get_statistic_value();
+	fit.set_param_value(pname,lower);
+	fit.fit();
+	Ts s=fit.get_statistic_value();
+	fit.set_param_value(pname,(lower+lower1)/2);
+	fit.fit();
+	Ts s01=fit.get_statistic_value();
+	//std::cerr<<"l:"<<lower1<<"\t"<<(lower1+lower)/2<<"\t"<<lower<<std::endl;
+	std::cerr<<s1-target_chi<<"\t"<<s01-target_chi<<"\t"<<s-target_chi<<std::endl;
+	if((s1-target_chi)*(s01-target_chi)<0&&(s-target_chi)*(s01-target_chi)>=0)
+	  {
+	    lower=(lower+lower1)/2;
+	  }
+	else if((s1-target_chi)*(s01-target_chi)>0&&(s-target_chi)*(s01-target_chi)<=0)
+	  {
+	    lower1=(lower+lower1)/2;
+	  }
+	else
+	  {
+	    std::cerr<<"strange statistic structure"<<std::endl;
+	    break;
+	  }
+      }
+
+    dynamic_cast<freeze_param<Ty,Tx,Tp,Tstr>& >(fit.get_param_modifier())-=freeze_param<Ty,Tx,Tp,Tstr>(pname);
+    
+  }
+}
+
+
+#endif
+//EOF
-- 
cgit v1.2.2