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 ++++++++++++++++++++++++++++++++++++++--------------
 misc/optvec.hpp      | 32 ++++++++++++++++++---
 statistics/chisq.hpp | 18 ++++++++++--
 test/test_fitter.cpp |  4 +++
 vmodels/beta1d.hpp   | 64 +++++++++++++++++++++++++++++++++++++++++
 vmodels/constant.hpp | 54 +++++++++++++++++++++++++++++++++++
 vmodels/poly1d.hpp   | 78 ++++++++++++++++++++++++++++++++++++++++++++++++++
 7 files changed, 302 insertions(+), 28 deletions(-)
 create mode 100644 vmodels/beta1d.hpp
 create mode 100644 vmodels/constant.hpp
 create mode 100644 vmodels/poly1d.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)),
diff --git a/misc/optvec.hpp b/misc/optvec.hpp
index cc85ccf..b259029 100644
--- a/misc/optvec.hpp
+++ b/misc/optvec.hpp
@@ -25,8 +25,8 @@ namespace opt_utilities
 
     optvec& operator=(const optvec& rhs)
     {
-      //dynamic_cast<std::vector<T>&>(*this).operator=(rhs);
-      std::vector<T>::operator=(rhs);
+      static_cast<std::vector<T>&>(*this).operator=(rhs);
+      //std::cerr<<"rhs.size="<<rhs.size()<<std::endl;
       return *this;
     }
   public:
@@ -208,8 +208,8 @@ namespace opt_utilities
   optvec<T> operator/(const T& x1,const optvec<T>& x2)
   {
     
-    optvec<T> result(x1.size());
-    for(size_t i=0;i!=x1.size();++i)
+    optvec<T> result(x2.size());
+    for(size_t i=0;i!=x2.size();++i)
       {
 	result[i]=x1/x2[i];
       }
@@ -260,6 +260,30 @@ namespace opt_utilities
     return false;
   }
 
+  template <typename Ty>
+  inline Ty randn(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;
+  }
+  
+  
+  template <typename T>
+  optvec<T> rand_norm(const optvec<T>& y0, const optvec<T>& y_err)
+  {
+    optvec<T> result(y0.size());
+    for(int i=0;i<y0.size();++i)
+      {
+	result[i]=randn(y0[i],y_err[i]);
+      }
+    return result;
+  }
+
 }
 
 
diff --git a/statistics/chisq.hpp b/statistics/chisq.hpp
index 339ab6f..c460394 100644
--- a/statistics/chisq.hpp
+++ b/statistics/chisq.hpp
@@ -247,14 +247,26 @@ namespace opt_utilities
       Ts result(0);
       for(int i=(this->get_data_set()).size()-1;i>=0;--i)
 	{
-	  Ty chi=(this->get_data_set().get_data(i).get_y()-eval_model(this->get_data_set().get_data(i).get_x(),p))/this->get_data_set().get_data(i).get_y_upper_err();
+	  Ty chi(this->get_data_set().get_data(0).get_y().size());
+	  for(int j=0;j<chi.size();++j)
+	    {
+	      Ty model_y(eval_model(this->get_data_set().get_data(i).get_x(),p));
+	      if(model_y[j]>this->get_data_set().get_data(i).get_y()[j])
+		{
+		  chi[j]=(this->get_data_set().get_data(i).get_y()[j]-model_y[j])/this->get_data_set().get_data(i).get_y_upper_err()[j];
+		}
+	      else
+		{
+		  chi[j]=(this->get_data_set().get_data(i).get_y()[j]-model_y[j])/this->get_data_set().get_data(i).get_y_lower_err()[j];
+		}
+	    }
 	  result+=sum(chi*chi);
-
+	  
 	}
       return result;
     }
   };
-
+  
   
 }
 
diff --git a/test/test_fitter.cpp b/test/test_fitter.cpp
index 5e5b139..226d4ce 100644
--- a/test/test_fitter.cpp
+++ b/test/test_fitter.cpp
@@ -1,4 +1,5 @@
 #include <core/fitter.hpp>
+#include <misc/bootstrap.hpp>
 #include <methods/powell/powell_method.hpp>
 #include <statistics/chisq.hpp>
 #include <models/nbeta1d.hpp>
@@ -36,4 +37,7 @@ int main(int argc,char* argv[])
     {
       cout<<f.get_param_info(i).get_name()<<"="<<f.get_param_info(i).get_value()<<endl;
     }
+  bootstrap<double,double,vector<double>,double,std::string> bst;
+  bst.set_fitter(f);
+  bst.sample(1000);
 }
diff --git a/vmodels/beta1d.hpp b/vmodels/beta1d.hpp
new file mode 100644
index 0000000..42be3fd
--- /dev/null
+++ b/vmodels/beta1d.hpp
@@ -0,0 +1,64 @@
+#ifndef VBETA_MODEL_H_
+#define VBETA_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <misc/optvec.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class beta1d
+    :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+  {
+    typedef optvec<T> Tv;
+  private:
+    beta1d<T>* do_clone()const
+    {
+      return new beta1d<T>(*this);
+    }
+
+    const char* do_get_type_name()const
+    {
+      return "1d surface brightness beta model";
+    }
+  public:
+    beta1d()
+    {
+      this->push_param_info(param_info<Tv>("S0",1));
+      this->push_param_info(param_info<Tv>("rc",10));
+      this->push_param_info(param_info<Tv>("beta",2./3.));
+      this->push_param_info(param_info<Tv>("bkg",0));
+    }
+
+  public:
+    Tv do_eval(const Tv& x,const Tv& param)
+    {
+      Tv result(x.size());
+      T S0=get_element(param,0);
+      T r_c=get_element(param,1);
+      T beta=get_element(param,2);
+      T bkg=get_element(param,3);
+      
+      //return x*get_element(param,0)+get_element(param,1);
+      for(size_t i=0;i!=x.size();++i)
+	{
+	  
+	  result[i]=bkg+S0*pow(1+(x[i]*x[i])/(r_c*r_c),-3*beta+static_cast<T>(.5));
+
+	}
+      return result;
+    }
+
+  private:
+    std::string do_get_information()const
+    {
+      return "";
+    }
+  };
+}
+
+
+
+#endif
+//EOF
diff --git a/vmodels/constant.hpp b/vmodels/constant.hpp
new file mode 100644
index 0000000..c4044b0
--- /dev/null
+++ b/vmodels/constant.hpp
@@ -0,0 +1,54 @@
+#ifndef VCONSTANT_MODEL_H_
+#define VCONSTANT_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <misc/optvec.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class constant
+    :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+  {
+    typedef optvec<T> Tv;
+  private:
+    constant<T>* do_clone()const
+    {
+      return new constant<T>(*this);
+    }
+    const char* do_get_type_name()const
+    {
+      return "constant";
+    }
+  public:
+    constant()
+    {
+      this->push_param_info(param_info<Tv>("c",1));
+    }
+
+  public:
+    Tv do_eval(const Tv& x,const Tv& param)
+    {
+      //return x*param[0]+param[1];
+      Tv result(x.size());
+      for(int i=0;i<result.size();++i)
+	{
+	  result[i]=param[0];
+	}
+      return result;
+    }
+
+  private:
+    std::string do_get_information()const
+    {
+      return "Constant\n"
+	"y=C\n";
+    }
+  };
+}
+
+
+
+#endif
+//EOF
diff --git a/vmodels/poly1d.hpp b/vmodels/poly1d.hpp
new file mode 100644
index 0000000..adaf03b
--- /dev/null
+++ b/vmodels/poly1d.hpp
@@ -0,0 +1,78 @@
+#ifndef POLY_MODEL_H_
+#define POLY_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <cmath>
+#include <sstream>
+#include <cassert>
+
+namespace opt_utilities
+{
+  template <typename T,int n>
+  class poly1d
+    :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+  {
+    typedef optvec<T> Tv;
+  private:
+    poly1d<T,n>* do_clone()const
+    {
+      return new poly1d<T,n>(*this);
+    }
+
+    const char* do_get_type_name()const
+    {
+      return "polynomial";
+    }
+  public:
+    poly1d()
+    {
+      assert(n>=0);
+      for(int i=0;i<=n;++i)
+	{
+	  std::ostringstream ostr;
+	  ostr<<"a"<<i;
+	  this->push_param_info(param_info<Tv>(ostr.str().c_str(),1));  
+	}
+      
+    }
+
+  public:
+    Tv do_eval(const Tv& x,const Tv& param)
+    {
+      //      return x*get_element(param,0)+get_element(param,1);
+      Tv result(x.size());
+      for(int m=0;m<result.size();++m)
+	{
+	  result[m]=0;
+	  for(int i=0;i<=n;++i)
+	    {
+	      typename Tv::value_type xn(1);
+	      for(int j=0;j<i;++j)
+		{
+		  xn*=x[m];
+		}
+	      result[m]+=get_element(param,i)*xn;
+	    }
+	}
+      return result;
+    }
+
+  private:
+    std::string do_get_information()const
+    {
+      std::ostringstream ostr;
+      ostr<<n<<"-order polynorminal model\n";
+      for(int i=0;i<n;++i)
+	{
+	  ostr<<"a"<<i<<"*x^"<<i<<"+";
+	}
+      ostr<<"a"<<n<<"*x^"<<n;
+      return ostr.str();
+    }
+  };
+}
+
+
+
+#endif
+//EOF
-- 
cgit v1.2.2