From 1f4a944064bc42284c33e6b755353d191cf288e8 Mon Sep 17 00:00:00 2001
From: astrojhgu <astrojhgu@ed2142bd-67ad-457f-ba7c-d818d4011675>
Date: Mon, 15 Dec 2008 07:26:12 +0000
Subject: git-svn-id: file:///home/svn/opt_utilities@1
 ed2142bd-67ad-457f-ba7c-d818d4011675

---
 core/default_data_set.hpp |   47 ++
 core/fitter.hpp           | 1043 +++++++++++++++++++++++++++++++++++++++++++++
 core/freeze_param.hpp     |  204 +++++++++
 core/opt_exception.hpp    |  106 +++++
 core/opt_traits.hpp       |   64 +++
 core/optimizer.hpp        |  289 +++++++++++++
 6 files changed, 1753 insertions(+)
 create mode 100644 core/default_data_set.hpp
 create mode 100644 core/fitter.hpp
 create mode 100644 core/freeze_param.hpp
 create mode 100644 core/opt_exception.hpp
 create mode 100644 core/opt_traits.hpp
 create mode 100644 core/optimizer.hpp

(limited to 'core')

diff --git a/core/default_data_set.hpp b/core/default_data_set.hpp
new file mode 100644
index 0000000..be06c04
--- /dev/null
+++ b/core/default_data_set.hpp
@@ -0,0 +1,47 @@
+#ifndef DEFAULT_DATA_SET
+#define DEFAULT_DATA_SET
+#include "fitter.hpp"
+#include <vector>
+
+
+namespace opt_utilities
+{
+
+template <typename Ty,typename Tx>
+class default_data_set
+  :public data_set<Ty,Tx>
+{
+private:
+  std::vector<data<Ty,Tx> > data_vec;
+
+  data_set<Ty,Tx>* do_clone()const
+  {
+    return new default_data_set<Ty,Tx>(*this);
+  }
+
+
+  const data<Ty,Tx>& do_get_data(int i)const
+  {
+    return data_vec.at(i);
+  }
+  
+  int do_size()const
+  {
+    return data_vec.size();
+  }
+  
+  void do_push_back(const data<Ty,Tx>& d)
+  {
+    data_vec.push_back(d);
+  }
+  
+  void do_clear()
+  {
+    data_vec.clear();
+  }
+  
+};
+}
+
+#endif
+//EOF
diff --git a/core/fitter.hpp b/core/fitter.hpp
new file mode 100644
index 0000000..ed90a83
--- /dev/null
+++ b/core/fitter.hpp
@@ -0,0 +1,1043 @@
+#ifndef FITTER_HPP
+#define FITTER_HPP
+#include "opt_exception.hpp"
+#include "optimizer.hpp"
+#include <vector>
+#include <string>
+#include <cstdlib>
+#include <cassert>
+#include <iostream>
+namespace opt_utilities
+{
+
+  ///////////////////////////////////
+  //////class data///////////////////
+  //////contain single data point////
+  ///////////////////////////////////
+  template<typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>
+  class statistic;
+
+  template<typename Ty,typename Tx,typename Tp,typename Tstr>
+  class param_modifier;
+
+
+  template<typename Ty,typename Tx>
+  class data
+  {
+  private:
+    Tx x,x_lower_err,x_upper_err;
+    Ty y,y_lower_err,y_upper_err;
+  public:
+    data(const Tx& _x,const Ty& _y,const Ty& _y_lower_err,const Ty& _y_upper_err,const Tx& _x_lower_err,const Tx& _x_upper_err)
+    {
+      opt_eq(x,_x);
+      opt_eq(x_lower_err,_x_lower_err);
+      opt_eq(x_upper_err,_x_upper_err);
+      opt_eq(y,_y);
+      opt_eq(y_lower_err,_y_lower_err);
+      opt_eq(y_upper_err,_y_upper_err);
+ 
+    }
+
+    data()
+      :x(),
+       x_lower_err(),
+       x_upper_err(),
+       y(),
+       y_lower_err(),
+       y_upper_err()
+    {}
+
+    data(const data& rhs)
+    {
+      opt_eq(x,rhs.x);
+      opt_eq(x_lower_err,rhs.x_lower_err);
+      opt_eq(x_upper_err,rhs.x_upper_err);
+      opt_eq(y,rhs.y);
+      opt_eq(y_lower_err,rhs.y_lower_err);
+      opt_eq(y_upper_err,rhs.y_upper_err);
+    }
+
+    data& operator=(const data& rhs)
+    {
+      opt_eq(x,rhs.x);
+      opt_eq(x_lower_err,rhs.x_lower_err);
+      opt_eq(x_upper_err,rhs.x_upper_err);
+      opt_eq(y,rhs.y);
+      opt_eq(y_lower_err,rhs.y_lower_err);
+      opt_eq(y_upper_err,rhs.y_upper_err);
+      return *this;
+    }
+
+  public:
+    const Tx& get_x()const
+    {
+      return x;
+    }
+    
+    const Tx& get_x_lower_err()const
+    {
+      return x_lower_err;
+    }
+    
+    const Tx& get_x_upper_err()const
+    {
+      return x_upper_err;
+    }
+
+    const Ty& get_y()const
+    {
+      return y;
+    }
+    
+    const Ty& get_y_lower_err()const
+    {
+      return y_lower_err;
+    }
+    
+    const Ty& get_y_upper_err()const
+    {
+      return y_upper_err;
+    }
+
+
+    void set_x(const Tx& _x)
+    {
+      opt_eq(x,_x);
+    }
+    
+    void set_x_lower_err(const Tx& _x)
+    {
+      opt_eq(x_lower_err,_x);
+    }
+    
+    void set_x_upper_err(const Tx& _x)
+    {
+      opt_eq(x_upper_err,_x);
+    }
+
+    void set_y(const Ty& _y)
+    {
+      opt_eq(y,_y);
+    }
+    
+    void set_y_lower_err(const Ty& _y)
+    {
+      opt_eq(y_lower_err,_y);
+    }
+    
+    void set_y_upper_err(const Ty& _y)
+    {
+      opt_eq(y_upper_err,_y);
+    }
+    
+
+  };
+
+  ////////////////////////////////////
+  ///class data_set///////////////////
+  ///contain a set of data////////////
+  ////////////////////////////////////
+  
+  
+  template <typename Ty,typename Tx>
+  class data_set
+  {
+  private:
+    
+    virtual const data<Ty,Tx>& do_get_data(int i)const=0;
+    virtual int do_size()const=0;
+    virtual void do_push_back(const data<Ty,Tx>&)=0;
+    virtual void do_clear()=0;
+    virtual data_set<Ty,Tx>* do_clone()const=0;
+
+    virtual void do_destroy()
+    {
+      delete this;
+    }
+  public:
+    const data<Ty,Tx>& get_data(int i)const
+    {
+      return this->do_get_data(i);
+    }
+    int size()const
+    {
+      return do_size();
+    }
+    void push_back(const data<Ty,Tx>& d)
+    {
+      return do_push_back(d);
+    }
+    void clear()
+    {
+      do_clear();
+    }
+    
+    data_set<Ty,Tx>* clone()const
+    {
+      return this->do_clone();
+    }
+
+    void destroy()
+    {
+      do_destroy();
+    }
+
+    virtual ~data_set(){}
+      
+  };
+
+  ///////////////////////////////////////////////
+  /////class param_info//////////////////////////
+  /////record the information of one parameter///
+  /////including the name, default value/////////
+  ///////////////////////////////////////////////
+
+  
+  template <typename Tp,typename Tstr=std::string>
+  class param_info
+  {
+  private:
+    Tstr name;
+    //bool frozen;
+    typename element_type_trait<Tp>::element_type default_value;
+    
+  public:
+    param_info(const Tstr& _name,
+	       const typename element_type_trait<Tp>::element_type& _v)
+      :name(_name),default_value(_v){}
+
+    param_info()
+      :name()
+    {}
+
+    param_info(const param_info& rhs)
+      :name(rhs.name)
+    {
+      opt_eq(default_value,rhs.default_value);
+    }
+    
+    param_info& operator=(const param_info& rhs)
+    {
+      name=rhs.name;
+      opt_eq(default_value,rhs.default_value);
+      return *this;
+    }
+
+    const Tstr& get_name()const
+    {
+      return this->name;
+    }
+
+    const typename element_type_trait<Tp>::element_type& get_default_value()const
+    {
+      return default_value;
+    }
+
+    void set_default_value(const typename element_type_trait<Tp>::element_type& x)
+    {
+      opt_eq(default_value,x);
+    }
+
+    void set_name(const Tstr& _name)
+    {
+      name=_name;
+    }
+  };
+
+
+
+
+
+  template <typename Ty,typename Tx,typename Tp,typename Tstr=std::string>
+  class model
+  {
+  private:
+    std::vector<param_info<Tp> > param_info_list;
+    param_info<Tp> null_param;
+    //    int num_free_params;
+    param_modifier<Ty,Tx,Tp,Tstr>* p_param_modifier;
+  private:
+    virtual model<Ty,Tx,Tp,Tstr>* do_clone()const=0;
+
+    virtual void do_destroy()
+    {
+      delete this;
+    }
+  public:
+    model<Ty,Tx,Tp,Tstr>* clone()const
+    {
+      return do_clone();
+    }
+
+    void destroy()
+    {
+      do_destroy();
+    }
+  public:
+    model()
+      :p_param_modifier(0)
+    {}
+    
+    model(const model& rhs)
+      :p_param_modifier(0)
+    {
+      if(rhs.p_param_modifier!=0)
+	{
+	  set_param_modifier(*(rhs.p_param_modifier));
+	}
+      param_info_list=rhs.param_info_list;
+      null_param=rhs.null_param;
+      
+    }
+    
+    model& operator=(const model& rhs)
+    {
+      if(this==&rhs)
+	{
+	  return *this;
+	}
+      if(rhs.p_param_modifier!=0)
+	{
+	  set_param_modifier(*(rhs.p_param_modifier));
+	}
+      param_info_list=rhs.param_info_list;
+      null_param=rhs.null_param;
+      return *this;
+    }
+    
+    virtual ~model()
+    {
+      if(p_param_modifier)
+	{
+	  //delete p_param_modifier;
+	  p_param_modifier->destroy();
+	}
+    }
+    
+    void set_param_modifier(const param_modifier<Ty,Tx,Tp,Tstr>& pm)
+    {
+      if(p_param_modifier!=0)
+	{
+	  //delete p_param_modifier;
+	  p_param_modifier->destroy();
+	}
+      p_param_modifier=pm.clone();
+      p_param_modifier->set_model(*this);
+    }
+
+    void set_param_modifier()
+    {
+      if(p_param_modifier!=0)
+	{
+	  //delete p_param_modifier;
+	  p_param_modifier->destroy();
+	}
+      p_param_modifier=0;
+    }
+    
+    param_modifier<Ty,Tx,Tp,Tstr>& get_param_modifier()
+    {
+      if(p_param_modifier==0)
+	{
+	  throw param_modifier_undefined();
+	}
+      return *p_param_modifier;
+    }
+
+    Tstr report_param_status(const Tstr& s)const
+    {
+      if(p_param_modifier==0)
+	{
+	  return Tstr();
+	}
+      
+      return p_param_modifier->report_param_status(s);
+	
+    }
+    
+
+    const param_info<Tp,Tstr>& get_param_info(const Tstr& pname)
+    {
+      for(typename std::vector<param_info<Tp,Tstr> >::iterator i=param_info_list.begin();
+	  i!=param_info_list.end();++i)
+	{
+	  if(i->get_name()==pname)
+	    {
+	      return *i;
+	    }
+	}
+      std::cerr<<"Param unfound!"<<std::endl;
+      assert(false);
+      throw param_not_found();
+      return null_param;
+    }
+
+    const param_info<Tp,Tstr>& get_param_info(int n)const
+    {
+      return param_info_list[n%get_num_params()];
+    }
+
+    
+    Tp get_all_params()const
+    {
+      Tp result;
+      resize(result,param_info_list.size());
+      for(size_t i=0;i<param_info_list.size();++i)
+	{
+	  //opt_eq(get_element(result,i),param_info_list[i].get_default_value());
+	  set_element(result,i,param_info_list[i].get_default_value());
+	  //get_element((Tp)result,i);
+	}
+      return result;
+    }
+
+  
+    int get_num_params()const
+    {
+      return (int)param_info_list.size();
+    }
+
+    int get_num_free_params()const
+    {
+      if(p_param_modifier)
+	{
+	  return (int)p_param_modifier->get_num_free_params();
+	}
+      return get_num_params();
+    }
+
+    void set_param_value(const Tstr& pname,
+			 const typename element_type_trait<Tp>::element_type& v)
+    {
+      //int porder=0;
+      for(typename std::vector<param_info<Tp,Tstr> >::iterator i=param_info_list.begin();
+	  i!=param_info_list.end();++i)
+	{
+	  if(i->get_name()==pname)
+	    {
+	      i->set_default_value(v);
+	      return;
+	    }
+	}
+      std::cerr<<"param "<<pname<<" unfound"<<std::endl;
+      throw param_not_found();
+    }
+
+    void set_param_value(const Tp& param)
+    {
+      for(size_t i=0;i<param_info_list.size();++i)
+	{
+	  param_info_list[i].set_default_value(get_element(param,i));
+	}
+    }
+
+    int get_param_order(const Tstr& pname)const
+    {
+      for(int i=0;i<(int)param_info_list.size();++i)
+	{
+	  if(param_info_list[i].get_name()==pname)
+	    {
+	      return i;
+	    }
+	}
+      //      std::cerr<<"param name invalid!"<<std::endl;
+      // assert(false);
+      throw param_not_found();
+      return -1;
+    }
+    
+
+		
+  protected:
+    void push_param_info(const param_info<Tp,Tstr>& pinfo)
+    {
+      param_info_list.push_back(pinfo);
+      //      this->num_free_params++;
+    }
+    
+    void clear_param_info()
+    {
+      //      this->num_free_params=0;
+      param_info_list.clear();
+    }
+
+
+
+  public:
+    Tstr to_string()const
+    {
+      return do_to_string();
+    }
+
+    Tp reform_param(const Tp& p)const
+    {
+      if(p_param_modifier==0)
+	{
+	  return p;
+	}
+      return p_param_modifier->reform(p);
+    }
+
+    Tp deform_param(const Tp& p)const
+    {
+      if(p_param_modifier==0)
+	{
+	  return p;
+	}
+      return p_param_modifier->deform(p);
+    }
+
+    
+    Ty eval(const Tx& x,const Tp& p)
+    {
+      return do_eval(x,reform_param(p));
+    }
+
+    virtual Ty do_eval(const Tx& x,const Tp& p)=0;
+
+  private:
+    virtual Tstr do_to_string()const
+    {
+      return Tstr();
+    }
+    
+  };
+
+
+
+  template<typename Ty,typename Tx,typename Tp,typename Ts=Ty,typename Tstr=std::string>
+  class fitter
+  {
+  public:
+    model<Ty,Tx,Tp,Tstr>* p_model;
+    statistic<Ty,Tx,Tp,Ts,Tstr>* p_statistic;
+    data_set<Ty,Tx>* p_data_set;
+    optimizer<Ts,Tp> optengine;
+  public:
+    fitter()
+      :p_model(0),p_statistic(0),p_data_set(0),optengine()
+    {}
+
+    
+    fitter(const fitter& rhs)
+      :p_model(0),p_statistic(0),p_data_set(0),optengine()
+    {
+      if(rhs.p_model!=0)
+	{
+	  set_model(*(rhs.p_model));
+	}
+      if(rhs.p_statistic!=0)
+	{
+	  set_statistic(*(rhs.p_statistic));
+	  assert(p_statistic->p_fitter!=0);
+	}
+      if(rhs.p_data_set!=0)
+	{
+	  load_data(*(rhs.p_data_set));
+	}
+      optengine=rhs.optengine;
+    }
+    
+    fitter& operator=(const fitter& rhs)
+    {
+      if(this==&rhs)
+	{
+	  return *this;
+	}
+      if(rhs.p_model!=0)
+	{
+	  set_model(*(rhs.p_model));
+	}
+      if(rhs.p_statistic!=0)
+	{
+	  set_statistic(*(rhs.p_statistic));
+	}
+      if(rhs.p_data_set!=0)
+	{
+	  load_data(*(rhs.p_data_set));
+	}
+      optengine=rhs.optengine;
+      return *this;
+    }
+    
+    virtual ~fitter()
+    {
+      if(p_model!=0)
+	{
+	  //delete p_model;
+	  p_model->destroy();
+	}
+      if(p_statistic!=0)
+	{
+	  //delete p_statistic;
+	  p_statistic->destroy();
+	}
+      if(p_data_set!=0)
+	{
+	  //delete p_data_set;
+	  p_data_set->destroy();
+	}
+    }
+    
+
+    Ty eval_model(const Tx& x,const Tp& p)
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      return p_model->eval(x,p);
+    }
+    
+
+  public:
+    void set_model(const model<Ty,Tx,Tp,Tstr>& m)
+    {
+      if(p_model!=0)
+	{
+	  //delete p_model;
+	  p_model->destroy();
+	}
+      p_model=m.clone();
+      //p_model=&m;
+      //  current_param.resize(m.get_num_params());
+    }
+    /*
+    void set(const model<Ty,Tx,Tp>& m)
+    {
+      set_model(m);
+    }
+    */
+    void set_statistic(const statistic<Ty,Tx,Tp,Ts,Tstr>& s)
+    {
+      if(p_statistic!=0)
+	{
+	  //delete p_statistic;
+	  p_statistic->destroy();
+	}
+      p_statistic=s.clone();
+      //p_statistic=&s;
+      p_statistic->set_fitter(*this);
+    }
+    /*
+    void set(const statistic<Ty,Tx,Tp>& s)
+    {
+      set_statistic(s);
+    }
+    */
+
+    void set_param_modifier(const param_modifier<Ty,Tx,Tp,Tstr>& pm)
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      p_model->set_param_modifier(pm);
+    }
+
+    void set_param_modifier()
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      p_model->set_param_modifier();
+    }
+
+    param_modifier<Ty,Tx,Tp,Tstr>& get_param_modifier()
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      return p_model->get_param_modifier();
+    }
+
+    Tstr report_param_status(const Tstr& s)const
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      return p_model->report_param_status(s);
+    }
+    
+    /*
+    void set(const param_modifier<Ty,Tx,Tp>& pm)
+    {
+      set_param_modifier(pm);
+    }
+    */
+    void load_data(const data_set<Ty,Tx>& da)
+    {
+      if(p_data_set!=0)
+	{
+	  //delete p_data_set;
+	  p_data_set->destroy();
+	}
+      p_data_set=da.clone();
+      if(p_statistic!=0)
+	{
+	  p_statistic->set_fitter(*this);
+	}
+    }
+	
+    const data_set<Ty,Tx>& datas()const
+    {
+      if(p_data_set==0)
+	{
+	  throw data_unloaded();
+	}
+      return *(this->p_data_set);
+    }
+
+    model<Ty,Tx,Tp,Tstr>& model()
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      return *(this->p_model);
+    }
+
+    statistic<Ty,Tx,Tp,Ts,Tstr>& statistic()
+    {
+      if(p_statistic==0)
+	{
+	  throw statistic_undefined();
+	}
+      return *(this->p_statistic);
+    }
+
+    opt_method<Ts,Tp>& method()
+    {
+      return optengine.method();
+    }
+    
+
+  public:
+    void set_param_value(const Tstr& pname,
+			 const typename element_type_trait<Tp>::element_type& v)
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      p_model->set_param_value(pname,v);
+    }
+
+    void set_param_value(const Tp& param)
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      p_model->set_param_value(param);
+    }
+
+    typename element_type_trait<Tp>::element_type get_param_value(const Tstr& pname)const
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      return p_model->get_param_info(pname).get_default_value();
+    }
+
+    const param_info<Tp,Tstr>& get_param_info(const Tstr& pname)const
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      return p_model->get_param_info(pname);
+    }
+
+    const param_info<Tp,Tstr>& get_param_info(int n)const
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      return p_model->get_param_info(n);
+    }
+
+    const int get_param_order(const Tstr& pname)const
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      return p_model->get_param_order(pname);
+    }
+
+    int get_num_params()const
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      return p_model->get_num_params();
+    }
+
+    Tp get_all_params()const
+    {
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      //return current_param;
+      return p_model->get_all_params();
+    }
+
+    void set_method(const opt_method<Ts,Tp>& pm)
+    {
+      //assert(p_optimizer!=0);
+      optengine.set_opt_method(pm);
+    }
+    /*
+    void set(const opt_method<Ty,Tp>& pm)
+    {
+      set_method(pm);
+    }
+    */
+    void set_precision(typename element_type_trait<Tp>::element_type y)
+    {
+      optengine.set_precision(y);
+    }
+
+    Tp fit()
+    {
+      //      assert(p_model!=0);
+      if(p_model==0)
+	{
+	  throw model_undefined();
+	}
+      if(p_data_set==0)
+	{
+	  throw data_unloaded();
+	}
+      //assert(p_optimizer!=0);
+      //assert(p_data_set!=0);
+      //assert(p_statistic!=0);
+      if(p_statistic==0)
+	{
+	  throw statistic_undefined();
+	}
+
+      optengine.set_func_obj(*p_statistic);
+      Tp current_param;
+      opt_eq(current_param,p_model->get_all_params());
+      Tp start_point;
+      opt_eq(start_point,p_model->deform_param(current_param));
+      //      std::cout<<start_point.size()<<std::endl;
+      
+      
+      //for(int i=0;i<(int)start_point.size();++i)
+      //	{
+      //	  std::cout<<start_point[i]<<",";
+      //	}
+      //std::cout<<std::endl;
+      //assert(start_point.size()!=0);
+      if(get_size(start_point)==0)
+	{
+	  //return start_point;
+	  return p_model->get_all_params();
+	}
+      optengine.set_start_point(start_point);
+      
+      Tp result;
+      opt_eq(result,optengine.optimize());
+
+      Tp decurrent_param;
+      opt_eq(decurrent_param,p_model->reform_param(result));
+      //current_param.resize(decurrent_param.size());
+      resize(current_param,get_size(decurrent_param));
+      opt_eq(current_param,decurrent_param);
+      p_model->set_param_value(current_param);
+      //   return current_param;
+      return p_model->get_all_params();
+    }
+
+  };
+
+
+  template<typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr=std::string>
+  class statistic
+    :public func_obj<Ts,Tp>
+  {
+  public:
+    fitter<Ty,Tx,Tp,Ts,Tstr>* p_fitter;
+  
+  private:
+    virtual statistic<Ty,Tx,Tp,Ts,Tstr>* do_clone()const=0;
+
+    virtual void do_destroy()
+    {
+      delete this;
+    }
+
+  public:
+    statistic<Ty,Tx,Tp,Ts,Tstr>* clone()const
+    {
+      return this->do_clone();
+    }
+    
+    void destroy()
+    {
+      return do_destroy();
+    }
+    statistic()
+      :p_fitter(0)
+    {}
+    
+    statistic(const statistic& rhs)
+      :p_fitter(rhs.p_fitter)
+    {}
+
+    statistic& operator=(const statistic& rhs)
+    {
+      if(this==&rhs)
+	{
+	  return *this;
+	}
+      p_fitter=rhs.p_fitter;
+      return *this;
+    }
+    
+    virtual ~statistic()
+    {}
+
+    virtual void set_fitter(fitter<Ty,Tx,Tp,Ts,Tstr>& pfitter)
+    {
+      p_fitter=&pfitter;
+    }
+
+    virtual const fitter<Ty,Tx,Tp,Ts,Tstr>& get_fitter()const
+    {
+      if(p_fitter==0)
+	{
+	  throw fitter_unset();
+	}
+      return *p_fitter;
+    }
+		
+    Ty eval_model(const Tx& x,const Tp& p)
+    {
+      if(p_fitter==0)
+	{
+	  throw fitter_unset();
+	}
+      return p_fitter->eval_model(x,p);
+    }
+
+		
+    const data_set<Ty,Tx>& datas()const
+    {
+      if(p_fitter==0)
+	{
+	  throw fitter_unset();
+	}
+      return p_fitter->datas();
+    }
+
+  };
+
+  template <typename Ty,typename Tx,typename Tp,typename Tstr=std::string>
+  class param_modifier
+  {
+  private:
+    model<Ty,Tx,Tp,Tstr>* p_model;
+  public:
+    Tp reform(const Tp& p)const
+    {
+      return do_reform(p);
+    }
+    Tp deform(const Tp& p)const
+    {
+      return do_deform(p);
+    }
+
+    param_modifier<Ty,Tx,Tp,Tstr>* clone()const
+    {
+      return do_clone();
+    }
+
+    void destroy()
+    {
+      do_destroy();
+    }
+
+  public:
+
+    param_modifier()
+      :p_model(0)
+    {}
+    
+    param_modifier(const param_modifier& rhs)
+      :p_model(rhs.p_model)
+    {}
+
+    param_modifier& operator=(const param_modifier& rhs)
+    {
+      if(this==&rhs)
+	{
+	  return *this;
+	}
+      p_model=rhs.p_model;
+      return *this;
+    }
+    
+  public:
+    void set_model(model<Ty,Tx,Tp,Tstr>& pf)
+    {
+      p_model=&pf;
+      update();
+    }
+
+    const model<Ty,Tx,Tp,Tstr>& get_model()const
+    {
+      if(p_model==0)
+	{
+	  std::cout<<"dajf;asdjfk;";
+	  throw model_undefined();
+	}
+      return *(this->p_model);
+    }
+
+    int get_num_free_params()const
+    {
+      return do_get_num_free_params();
+    }
+
+    Tstr report_param_status(const Tstr& name)const
+    {
+      return do_report_param_status(name);
+    }
+
+
+    virtual ~param_modifier(){}
+  private:
+    virtual Tp do_reform(const Tp& p)const=0;
+    virtual Tp do_deform(const Tp& p)const=0;
+    virtual int do_get_num_free_params()const=0;
+    virtual Tstr do_report_param_status(const Tstr&)const=0;
+    virtual void update(){}
+
+    virtual param_modifier<Ty,Tx,Tp,Tstr>* do_clone()const=0;
+    
+    virtual void do_destroy()
+    {
+      delete this;
+    }
+
+  };
+  
+
+  
+};
+
+
+#endif
+//EOF
diff --git a/core/freeze_param.hpp b/core/freeze_param.hpp
new file mode 100644
index 0000000..8ba677a
--- /dev/null
+++ b/core/freeze_param.hpp
@@ -0,0 +1,204 @@
+#ifndef FREEZE_PARAM_HPP
+#define FREEZE_PARAM_HPP
+#include "fitter.hpp"
+#include <vector>
+#include <set>
+
+namespace opt_utilities
+{
+  template <typename Ty,typename Tx,typename Tp,typename Tstr=std::string>
+  class freeze_param
+    :public param_modifier<Ty,Tx,Tp,Tstr>
+  {
+  private:
+    std::set<Tstr> param_names;
+    std::vector<int> param_num;
+    int num_free;
+    
+  public:
+    freeze_param()
+    {
+      
+    }
+
+    freeze_param(const Tstr& name)
+    {
+      param_names.insert(name);
+    }
+    
+  private:
+    freeze_param<Ty,Tx,Tp,Tstr>* do_clone()const
+    {
+      return new freeze_param<Ty,Tx,Tp,Tstr>(*this);
+    }
+    
+    
+
+    void update()
+    {
+      param_num.clear();
+      for(typename std::set<Tstr>::const_iterator  i=param_names.begin();
+	  i!=param_names.end();++i)
+	{
+	  try
+	    {
+	      param_num.push_back(this->get_model().get_param_order(*i));
+	    }
+	  catch(opt_exception& e)
+	    {
+	      param_names.erase(*i);
+	      throw;
+	    }
+	  
+	}
+    }
+
+    int do_get_num_free_params()const
+    {
+      return this->get_model().get_num_params()-param_num.size();
+    }
+
+    bool is_frozen(int i)const
+    {
+      if(find(param_num.begin(),param_num.end(),i)==param_num.end())
+	{
+	  return false;
+	}
+      return true;
+    }
+
+
+    Tp do_reform(const Tp& p)const
+    {
+      int nparams=(this->get_model().get_num_params());
+      Tp reformed_p(nparams);
+      int i=0,j=0;
+      for(i=0;i<(int)nparams;++i)
+	{
+	  if(this->is_frozen(i))
+	    {
+	      const param_info<Tp,Tstr>& pinf=this->get_model().get_param_info(i);
+	      //std::cout<<"frozen:"<<pinf.get_name()
+	      //	       <<i<<"\t"<<j
+	      //       <<std::endl;
+	      //opt_eq(get_element(reformed_p,i),pinf.get_default_value());
+	      set_element(reformed_p,i,pinf.get_default_value());
+	    }
+	  else
+	    {
+	      //opt_eq(get_element(reformed_p,i),get_element(p,j));
+	      set_element(reformed_p,i,get_element(p,j));
+	      j++;
+	    }
+	}
+      /*
+      for(int i=0;i<reformed_p.size();++i)
+      {
+      std::cout<<get_element(reformed_p,i)<<",";
+      }
+      */
+      //std::cout<<"\n";
+      return reformed_p;
+      // return p;
+    }
+
+    Tp do_deform(const Tp& p)const
+    {
+      Tp deformed_p(do_get_num_free_params());
+      int i(0),j(0);
+      assert(get_size(p)==this->get_model().get_num_params());
+      for(;i<(int)get_size(p);++i)
+	{
+	  //std::cout<<is_frozen(j)<<"\n";
+	  if(!this->is_frozen(i))
+	    {
+	      //opt_eq(get_element(deformed_p,j),get_element(p,i));
+	      set_element(deformed_p,j,get_element(p,i));
+	      j++;
+	    }
+	}
+      
+      assert(j==do_get_num_free_params());
+      return deformed_p;
+    }
+
+
+    Tstr do_report_param_status(const Tstr& name)const
+    {
+      if(param_names.find(name)==param_names.end())
+	{
+	  return "thawed";
+	}
+      return "frozen";
+    }
+
+  public:
+    freeze_param operator+(const freeze_param& fp)const
+    {
+      freeze_param result(*this);
+      for(typename std::set<Tstr>::const_iterator i=fp.param_names.begin();
+	  i!=fp.param_names.end();
+	  ++i)
+	{
+	  result.param_names.insert(*i);
+	}
+      return result;
+    }
+
+    freeze_param& operator+=(const freeze_param& fp)
+    {
+      //param_names.insert(param_names.end(),
+      //fp.param_names.begin(),
+      //fp.param_names.end());
+      for(typename std::set<Tstr>::const_iterator i=fp.param_names.begin();
+	  i!=fp.param_names.end();
+	  ++i)
+	{
+	  param_names.insert(*i);
+	}
+      try
+	{
+	  update();
+	}
+      catch(opt_exception& e)
+	{
+	  throw;
+	}
+      return *this;
+    }
+
+    freeze_param& operator-=(const freeze_param& fp)
+    {
+      //param_names.insert(param_names.end(),
+      //fp.param_names.begin(),
+      //fp.param_names.end());
+      for(typename std::set<Tstr>::const_iterator i=fp.param_names.begin();
+	  i!=fp.param_names.end();
+	  ++i)
+	{
+	  param_names.erase(*i);
+	}
+      try
+	{
+	  update();
+	}
+      catch(opt_exception& e)
+	{
+	  throw;
+	}
+      return *this;
+    }
+
+  };
+
+  template <typename Ty,typename Tx,typename Tp,typename Tstr>
+  freeze_param<Ty,Tx,Tp,Tstr> freeze(const Tstr& name)
+  {
+    return freeze_param<Ty,Tx,Tp,Tstr>(name);
+  }
+
+};
+
+
+#endif
+//EOF
diff --git a/core/opt_exception.hpp b/core/opt_exception.hpp
new file mode 100644
index 0000000..925e4b5
--- /dev/null
+++ b/core/opt_exception.hpp
@@ -0,0 +1,106 @@
+#ifndef OPT_EXCEPTION
+#define OPT_EXCEPTION
+#include <string>
+#include <exception>
+namespace opt_utilities
+{
+  class opt_exception
+    :public std::exception
+  {
+  private:
+    std::string _what;
+  public:
+    opt_exception()
+    {};
+    
+    ~opt_exception()throw()
+    {}
+    
+    opt_exception(const std::string& str)
+      :_what(str)
+    {}
+
+    const char* what()const throw()
+    {
+      return _what.c_str();
+    }
+  };
+  
+  class target_function_undefined
+    :public opt_exception
+  {
+  public:
+    target_function_undefined()
+      :opt_exception("target function undefined")
+    {}
+  };
+
+  class opt_method_undefined
+    :public opt_exception
+  {
+  public:
+    opt_method_undefined()
+      :opt_exception("opt method undefined")
+    {}
+  };
+  
+  class fitter_unset
+    :public opt_exception
+  {
+  public:
+    fitter_unset()
+      :opt_exception("fitter_unset")
+    {}
+  };
+
+  class model_undefined
+    :public opt_exception
+  {
+  public:
+    model_undefined()
+      :opt_exception("model_undefined")
+    {}
+  };
+
+  class data_unloaded
+    :public opt_exception
+  {
+  public:
+    data_unloaded()
+      :opt_exception("data not loaded")
+    {}
+  };
+
+
+  class statistic_undefined
+    :public opt_exception
+  {
+  public:
+    statistic_undefined()
+      :opt_exception("statistic undefined")
+    {}
+  };
+  
+  class param_not_found
+    :public opt_exception
+  {
+  public:
+    param_not_found()
+      :opt_exception("param name invalid")
+    {}
+  };
+  
+  class param_modifier_undefined
+    :public opt_exception
+  {
+  public:
+    param_modifier_undefined()
+      :opt_exception("param modifier undefined")
+    {}
+  };
+
+};
+
+
+#endif
+//EOF
diff --git a/core/opt_traits.hpp b/core/opt_traits.hpp
new file mode 100644
index 0000000..c539a63
--- /dev/null
+++ b/core/opt_traits.hpp
@@ -0,0 +1,64 @@
+#ifndef ARRAY_OPERATION
+#define ARRAY_OPERATION
+#include <cstddef>
+namespace opt_utilities
+{
+  /////////Useful function///////////////////////////////////
+  template <typename T>
+  inline size_t get_size(const T& x)
+  {
+    return x.size();
+  }
+  
+  template <typename T>
+  class element_type_trait
+  {
+  public:
+    typedef typename T::value_type element_type;
+  };
+
+  template <typename T>
+  class return_type_trait
+  {
+  public:
+    typedef T value_type;
+    typedef T& reference_type;
+    typedef const T& const_reference_type;
+  };
+  
+  template <typename T>
+  inline typename return_type_trait<typename element_type_trait<T>::element_type>::const_reference_type get_element(const T& x,size_t i)
+  {
+    return x[i];
+  }
+  /*
+  template <typename T>
+  inline typename element_type_trait<T>::element_type& get_element(T& x,size_t i)
+  {
+    return x[i];
+  }
+  */
+
+  template<typename T,typename TX>
+  inline void set_element(T& x,size_t i,
+			  const TX& v)
+  {
+    x[i]=v;
+  }
+
+  template <typename T>
+  inline void resize(T& x,size_t s)
+  {
+    x.resize(s);
+  }
+
+  template <typename Tl,typename Tr>
+  inline Tl& opt_eq(Tl& lhs,const Tr& rhs)
+  {
+    return (lhs=rhs);
+  }
+};
+
+
+
+#endif
diff --git a/core/optimizer.hpp b/core/optimizer.hpp
new file mode 100644
index 0000000..ca73775
--- /dev/null
+++ b/core/optimizer.hpp
@@ -0,0 +1,289 @@
+#ifndef OPTIMZER_H_
+#define OPTIMZER_H_
+//#define DEBUG
+#include <cstddef>
+#include "opt_traits.hpp"
+#include "opt_exception.hpp"
+#include <cstdlib>
+#include <functional>
+#ifdef DEBUG
+#include <iostream>
+using namespace std;
+#endif
+
+namespace opt_utilities
+{
+  /////////Forward declare///////////////////////////////////
+  template <typename rT,typename pT>
+  class optimizer;
+  
+  template <typename rT,typename pT>
+  class func_obj;
+  
+  template <typename rT,typename pT>
+  class opt_method;
+
+
+  //////////////Target Function/////////////////////
+  ///An eval function should be implemented/////////
+  ///The eval function return the function value////
+  ///which is wrapped by the func_obj///////////////
+  //////////////////////////////////////////////////
+  template <typename rT,typename pT>
+  class func_obj
+    :public std::unary_function<pT,rT>
+  {
+  private:
+    virtual rT do_eval(const pT&)=0;
+    virtual func_obj<rT,pT>* do_clone()const=0;
+    virtual void do_destroy()
+    {
+      delete this;
+    }
+
+  public:
+  public:
+    func_obj<rT,pT>* clone()const
+    {
+      return do_clone();
+    }
+
+    void destroy()
+    {
+      do_destroy();
+    }
+
+    rT operator()(const pT& p)
+    {
+      return do_eval(p);
+    }
+    
+
+    rT eval(const pT& p)
+    {
+      return do_eval(p);
+    };
+    virtual ~func_obj(){};
+    //    virtual XT walk(XT,YT)=0;
+  };
+  
+
+  ///////////////Optimization method//////////////////////
+  
+  template <typename rT,typename pT>
+  class opt_method
+  {
+  public:
+    virtual void do_set_optimizer(optimizer<rT,pT>&)=0;
+    virtual void do_set_precision(rT)=0;
+    virtual pT do_optimize()=0;
+    virtual void do_set_start_point(const pT& p)=0;
+    virtual opt_method<rT,pT>* do_clone()const=0;
+
+    virtual void do_destroy()
+    {
+      delete this;
+    }
+  public:
+    void  set_optimizer(optimizer<rT,pT>& op)
+      {
+	do_set_optimizer(op);
+      };
+
+    void set_precision(rT x)
+    {
+      do_set_precision(x);
+    }
+
+    void set_start_point(const pT& p)
+    {
+      do_set_start_point(p);
+    }
+    
+    pT optimize()
+      {
+	return do_optimize();
+      };
+    
+    opt_method<rT,pT>* clone()const
+    {
+      return do_clone();
+    }
+
+    void destroy()
+    {
+      do_destroy();
+    }
+
+    virtual ~opt_method(){};
+  };
+  
+  
+  ///////////Optimizer////////////////////////////////////
+  template <typename rT,typename pT>
+  class optimizer
+  {
+  public:
+
+  private:
+
+    ////////////pointer to an optimization method objection////////////
+    ////////////The optimization method implements a certain method ///
+    ////////////Currently only Mont-carlo method is implemented////////
+    opt_method<rT,pT>* p_opt_method;
+    func_obj<rT,pT>* p_func_obj;
+    
+  public:
+     optimizer()
+      :p_opt_method(0),p_func_obj(0)
+    {}
+
+    optimizer(func_obj<rT,pT>& fc,const opt_method<rT,pT>& om)
+      :p_func_obj(fc.clone()),p_opt_method(om.clone())
+    {
+      p_opt_method->set_optimizer(*this);
+    }
+    
+    optimizer(const optimizer& rhs)
+      :p_opt_method(0),p_func_obj(0)
+    {
+      if(rhs.p_func_obj!=0)
+	{
+	  set_func_obj(*(rhs.p_func_obj));
+	}
+      if(rhs.p_opt_method!=0)
+	{
+	  set_opt_method(*(rhs.p_opt_method));
+	}
+    }
+
+    optimizer& operator=(const optimizer& rhs)
+    {
+      if(this==&rhs)
+	{
+	  return *this;
+	}
+      if(rhs.p_func_obj!=0)
+	{
+	  set_func_obj(*(rhs.p_func_obj));
+	}
+      if(rhs.p_opt_method!=0)
+	{
+	  set_opt_method(*(rhs.p_opt_method));
+	}
+      return *this;
+    }
+    
+
+    virtual ~optimizer()
+    {
+      if(p_func_obj!=0)
+	{
+	  //delete p_func_obj;
+	  p_func_obj->destroy();
+	}
+      if(p_opt_method!=0)
+	{
+	  //delete p_opt_method;
+	  p_opt_method->destroy();
+	}
+    };
+
+  public:
+    ////////////Re-set target function object///////////////////////////
+    void set_func_obj(const func_obj<rT,pT>& fc)
+    {
+      if(p_func_obj!=0)
+	{
+	  //delete p_func_obj;
+	  p_func_obj->destroy();
+	}
+      p_func_obj=fc.clone();
+      if(p_opt_method!=0)
+	{
+	  p_opt_method->set_optimizer(*this);
+	}
+    }
+
+    ////////////Re-set optimization method//////////////////////////////
+    void set_opt_method(const opt_method<rT,pT>& om)
+    {
+      if(p_opt_method!=0)
+	{
+	  //delete p_opt_method;
+	  p_opt_method->destroy();
+	}
+      
+      p_opt_method=om.clone();
+      p_opt_method->set_optimizer(*this);
+    }
+
+    opt_method<rT,pT>& method()
+    {
+      if(p_opt_method==0)
+	{
+	  throw opt_method_undefined();
+	}
+      return *(this->p_opt_method);
+    }
+
+    void set_precision(rT x)
+    {
+      if(p_opt_method==0)
+	{
+	  throw opt_method_undefined();
+	}
+      p_opt_method->set_precision(x);
+    }
+
+    void set_start_point(const pT& x)
+    {
+      if(p_opt_method==0)
+	{
+	  throw opt_method_undefined();
+	}
+      p_opt_method->set_start_point(x);
+    }
+    
+    ////////////Just call the eval function in the target function object///
+    ////////////In case the pointer to a target function is uninitialed/////
+    ////////////a zero-value is returned////////////////////////////////////
+    rT eval(const pT& x)
+    {
+      if(p_func_obj==0)
+	{
+	  throw target_function_undefined();
+	}
+      return p_func_obj->eval(x);
+    }
+
+	
+
+    ////////////Just call the optimize() function in the optimization method//
+    ////////////If no optimization method is given, an zero-value is returned/
+    pT optimize()
+    {
+      if(p_opt_method==0)
+	{
+	  throw opt_method_undefined();
+	}
+      if(p_func_obj==0)
+	{
+	  throw target_function_undefined();
+	}
+      return p_opt_method->optimize();
+    }
+    
+    ////////////Function that offers the access to the target function object///
+    func_obj<rT,pT>* ptr_func_obj()
+    {
+      return p_func_obj;
+    }
+    
+  };
+};
+  
+#endif
+//EOF
+
+
-- 
cgit v1.2.2