From db1d374b7a24b8f3f74e17d706704d93f1e66d90 Mon Sep 17 00:00:00 2001
From: astrojhgu <astrojhgu@ed2142bd-67ad-457f-ba7c-d818d4011675>
Date: Mon, 8 Feb 2010 09:23:14 +0000
Subject: git-svn-id: file:///home/svn/opt_utilities@110
 ed2142bd-67ad-457f-ba7c-d818d4011675

---
 misc/bootstrap.hpp | 80 ++++++++++++++++++++++++++++++++++++++++--------------
 1 file changed, 59 insertions(+), 21 deletions(-)

(limited to 'misc/bootstrap.hpp')

diff --git a/misc/bootstrap.hpp b/misc/bootstrap.hpp
index f2ba525..88332ed 100644
--- a/misc/bootstrap.hpp
+++ b/misc/bootstrap.hpp
@@ -3,11 +3,12 @@
 */
 
 
-#ifndef BOOT_STRIP
-#define BOOT_STRIP
+#ifndef BOOT_STRAP
+#define BOOT_STRAP
 #define OPT_HEADER
 #include <core/fitter.hpp>
 #include <vector>
+#include <list>
 #include <cstdlib>
 #include <iostream>
 #include <utility>
@@ -16,6 +17,18 @@
 using std::cerr;
 namespace opt_utilities
 {
+  template <typename Ty>
+  inline Ty rand_norm(Ty y0,Ty y_err)
+  {
+    Ty y;
+    do
+      {
+	y=(rand()/(Ty)RAND_MAX-(Ty).5)*(10*y_err)+y0;
+      }
+    while(rand()/(Ty)RAND_MAX>exp(-(y-y0)*(y-y0)/(y_err*y_err)));
+    return y;
+  }
+
 
   /**
      \brief using bootstrap method to estimate confidence interval
@@ -29,23 +42,16 @@ namespace opt_utilities
   class bootstrap
   {
   private:
-    Ty rand_norm(Ty y0,Ty y_err)const
-    {
-      Ty y;
-      do
-	{
-	  y=(rand()/(Ty)RAND_MAX-(Ty).5)*(10*y_err)+y0;
-	}
-      while(rand()/(Ty)RAND_MAX>exp(-(y-y0)*(y-y0)/(y_err*y_err)));
-      return y;
-      
-    }
   public:
-    std::vector<Tp> param_pool;
+    typedef typename std::list<Tp>::const_iterator param_iterator;
+  public:
+    std::list<Tp> param_pool;
     default_data_set<Ty,Tx> current_data_set;
     default_data_set<Ty,Tx> origin_data_set;
     fitter<Ty,Tx,Tp,Ts,Tstr>* p_fitter;
     Tp origin_param;
+    mutable bool bstop;
+    
   public:
 
     /**
@@ -83,11 +89,12 @@ namespace opt_utilities
      */
     void reset()
     {
-      if(p_fitter!=NULL)
+      if(p_fitter!=0)
 	{
 	  p_fitter->load_data(origin_data_set);
 	  p_fitter->set_param_value(origin_param);
 	}
+      p_fitter=0;
     }
 
 
@@ -97,9 +104,10 @@ namespace opt_utilities
      */
     void sample(int n)
     {
+      bstop=false;
       if(p_fitter!=NULL)
 	{
-	  for(int i=0;i<n;++i)
+	  for(int i=0;i<n&&!bstop;++i)
 	    {
 	      sample();
 	    }
@@ -113,15 +121,45 @@ namespace opt_utilities
 	}
     }
 
+    void stop()
+    {
+      bstop=true;
+    }
+
     /**
        get the param of i-th sampling
        \param i the order of parameter
        \return the parameter
      */
-    const Tp& get_param(int i)const
+    const Tp& get_param(int n)const
     {
-      return param_pool.at(i);
+      //return param_pool.at(i);
+      int cnt=0;
+      for(typename std::list<Tp>::const_iterator i=param_pool.begin();
+	  i!=param_pool.end();++i)
+	{
+	  if(cnt==n)
+	    {
+	      return *i;
+	    }
+	  ++cnt;
+	}
+      throw opt_exception("excesses param_pool size");
+      return *param_pool.begin();
+    }
+
+    typename std::list<Tp>::const_iterator 
+    param_begin()const
+    {
+      return param_pool.begin();
+    }
+
+    typename std::list<Tp>::const_iterator 
+    param_end()const
+    {
+      return param_pool.end();
     }
+    
 
     int get_param_pool_size()const
     {
@@ -140,7 +178,7 @@ namespace opt_utilities
 	{
 	  data<Ty,Tx> d;
 	  d=origin_data_set.get_data(i);
-	  d.set_y(rand_norm(d.get_y(),(d.get_y_upper_err()+d.get_y_lower_err())/2));
+	  d.set_y(rand_norm(d.get_y(),(d.get_y_upper_err()+d.get_y_lower_err())/2.));
 	  current_data_set.add_data(d);
 	}
       p_fitter->load_data(current_data_set);
@@ -176,7 +214,7 @@ namespace opt_utilities
       //sample();
       std::vector<typename element_type_trait<Tp>::element_type> _tmp;
       int order=p_fitter->get_param_order(param_name);
-      for(typename std::vector<Tp>::iterator i=param_pool.begin();
+      for(typename std::list<Tp>::iterator i=param_pool.begin();
 	  i!=param_pool.end();++i)
 	{
 	  _tmp.push_back((*i)[order]);
@@ -186,7 +224,7 @@ namespace opt_utilities
 	typename std::vector<typename element_type_trait<Tp>::element_type>::iterator> 
 	itv=equal_range(_tmp.begin(),_tmp.end(),origin_param[order]);
       int current_param_position=itv.second-_tmp.begin();
-      std::cerr<<_tmp.size()<<std::endl;
+      //std::cerr<<_tmp.size()<<std::endl;
       return std::pair<typename element_type_trait<Tp>::element_type,
 	typename element_type_trait<Tp>::element_type>(
 				 _tmp.at((int)((1-level)*current_param_position)),
-- 
cgit v1.2.2