git-svn-id: file:///home/svn/opt_utilities@138 ed2142bd-67ad-457f-ba7c-d818d4011675

3 files changed, 312 insertions, 0 deletions

diff --git a/distributions/component.hpp b/distributions/component.hpp
new file mode 100644
index 0000000..b595cab
--- /dev/null
+++ b/distributions/component.hpp
@@ -0,0 +1,179 @@
+#ifndef COMPONENT_MODEL_H_
+#define COMPONENT_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <cmath>
+#include <misc/optvec.hpp>
+#include <sstream>
+#include <iostream>
+/*
+ *Represents a distribution composed of more than one components
+ */
+
+
+namespace opt_utilities
+{
+  template <typename T>
+  class component
+    :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+  {
+  private:
+    std::vector<model<optvec<T>,optvec<T>,optvec<T>,std::string>*> components;
+    std::vector<int> weight_num;
+  private:
+    component* do_clone()const
+    {
+      return new component<T>(*this);
+    }
+
+    const char* do_get_type_name()const
+    {
+      return "Multi-distribution";
+    }
+  public:
+    component()
+    {
+      //this->push_param_info(param_info<optvec<T> >("x0",0));
+      //this->push_param_info(param_info<optvec<T> >("sigma",1));
+    }
+
+    component(const component& rhs)
+    {
+      for(int i=0;i<rhs.components.size();++i)
+	{
+	  add_component(*rhs.components[i],rhs.get_param_info(rhs.weight_num[i]).get_value());
+	}
+      for(int i=0;i<rhs.get_num_params();++i)
+	{
+	  set_param_info(rhs.get_param_info(i));
+	}
+      
+    }
+
+    component& operator=(const component& rhs)
+    {
+      if(this==&rhs)
+	{
+	  return rhs;
+	}
+      for(int i=0;i<components.size();++i)
+	{
+	  components[i]->destroy();
+	}
+      components.clear();
+      for(int i=0;i<rhs.components.size();++i)
+	{
+	  add_component(*rhs.components[i]);
+	}
+      for(int i=0;i<rhs.get_num_params();++i)
+	{
+	  set_param_info(rhs.get_param_info(i));
+	}
+      
+    }
+
+    ~component()
+    {
+      for(int i=0;i<components.size();++i)
+	{
+	  components[i]->destroy();
+	}
+    }
+
+
+  public:
+    void add_component(const model<optvec<T>,optvec<T>,optvec<T>,std::string>& m,const T& w=0)
+    {
+      int morder=components.size();
+      components.push_back(m.clone());
+      std::ostringstream oss;
+      oss<<morder;
+      std::string smorder=oss.str();
+      if(components.size()>1)
+	{
+	  param_info<optvec<T>,std::string> p1(std::string("_w")+smorder,w);	  
+	  this->push_param_info(p1);
+	  //std::cout<<this->get_num_params()<<std::endl;
+	  weight_num.push_back(this->get_num_params()-1);
+	}
+
+      int np=m.get_num_params();
+      for(int i=0;i<np;++i)
+	{
+	  param_info<optvec<T>,std::string> p(m.get_param_info(i));
+	  param_info<optvec<T>,std::string> p1(p.get_name()+smorder,p.get_value());
+	  this->push_param_info(p1);
+	  if(this->get_num_params()==1)
+	    {
+	      std::cout<<"origin:"<<p.get_value()<<std::endl;
+	      std::cout<<"self:"<<this->get_param_info(0).get_value()<<std::endl;
+
+	    }
+	  //std::cout<<this->get_num_params()<<std::endl;
+	}
+    }
+
+    optvec<T> convert_unit_sph(const optvec<T>& p)
+    {
+      int ndim=p.size()+1;
+      optvec<double> result(ndim);
+      result[0]=1;
+      for(int i=0;i<p.size();++i)
+	{
+	  for(int j=0;j<=i;++j)
+	    {
+	      result[j]*=cos(p[i]);
+	    }
+	  result[i+1]=sin(p[i]);
+	}
+      for(int i=0;i<result.size();++i)
+	{
+	  result[i]=result[i]*result[i];
+	  //	  std::cout<<result[i]<<"\t";
+	}
+      //      std::cout<<std::endl;
+      return result;
+    }
+
+    
+    optvec<T> do_eval(const optvec<T>& x,const optvec<T>& param)
+    {
+      T result(0);
+      optvec<T> weight_angle;
+      for(int i=0;i<weight_num.size();++i)
+	{
+	  weight_angle.push_back(param[weight_num[i]]);
+	}
+      optvec<T> weight(convert_unit_sph(weight_angle));
+      
+      int pnum=0;
+      //std::cout<<param[2]<<"\t";
+      //std::cout<<weight[0]<<"\t"<<weight[1]<<std::endl;
+      for(int i=0;i<components.size();++i)
+	{
+	  T temp_result=0;
+	  optvec<T> p(components[i]->get_num_params());
+	  for(int j=0;j<p.size();++j)
+	    {
+	      p[j]=param[pnum++];
+	    }
+	  result+=(components[i]->eval(x,p)[0]*weight[i]);
+	  ++pnum;
+	}
+      optvec<T> result1(1);
+      result1[0]=result;;
+      return result1;
+    }
+
+  private:
+    std::string do_get_information()const
+    {
+      return "";
+    }
+  };
+}
+
+
+
+#endif
+//EOF
diff --git a/distributions/normed_dgauss1d.hpp b/distributions/normed_dgauss1d.hpp
new file mode 100644
index 0000000..4060fa6
--- /dev/null
+++ b/distributions/normed_dgauss1d.hpp
@@ -0,0 +1,75 @@
+#ifndef NDGAUSS_MODEL_H_
+#define NDGAUSS_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <cmath>
+#include <misc/optvec.hpp>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class normed_dgauss1d
+    :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+  {
+  private:
+    normed_dgauss1d* do_clone()const
+    {
+      return new normed_dgauss1d<T>(*this);
+    }
+
+    const char* do_get_type_name()const
+    {
+      return "1d double normed gaussian";
+    }
+  public:
+    normed_dgauss1d()
+    {
+      this->push_param_info(param_info<optvec<T> >("x01",0));
+      this->push_param_info(param_info<optvec<T> >("sigma1",1));
+      this->push_param_info(param_info<optvec<T> >("x02",0.1));
+      this->push_param_info(param_info<optvec<T> >("sigma2",1));
+      this->push_param_info(param_info<optvec<T> >("theta",1));      
+    }
+
+
+  public:
+    optvec<T> do_eval(const optvec<T>& x,const optvec<T>& param)
+    {
+      const double pi=3.14159265358979323846;
+      T x01=get_element(param,0);
+      T sigma1=get_element(param,1);
+      T x02=get_element(param,2);
+      T sigma2=get_element(param,3);
+      T theta=get_element(param,4);
+      if(sigma1*sigma1<std::numeric_limits<double>::epsilon())
+	{
+	  sigma1=std::numeric_limits<double>::epsilon();
+	}
+      if(sigma2*sigma2<std::numeric_limits<double>::epsilon())
+	{
+	  sigma2=std::numeric_limits<double>::epsilon();
+	}
+      T N1=1/sqrt(sigma1*sigma1*pi*2);
+      T N2=1/sqrt(sigma2*sigma2*pi*2);
+      
+      optvec<T> y1=(x-x01)/sigma1;
+      optvec<T> y2=(x-x02)/sigma2;
+      
+      T r1=sin(theta);
+      T r2=cos(theta);
+      
+      return r1*r1*N1*exp(-y1*y1/2.)+r2*r2*N2*exp(-y2*y2/2.);;
+    }
+
+  private:
+    std::string do_get_information()const
+    {
+      return "";
+    }
+  };
+}
+
+
+
+#endif
+//EOF
diff --git a/distributions/normed_gauss1d.hpp b/distributions/normed_gauss1d.hpp
new file mode 100644
index 0000000..09f5494
--- /dev/null
+++ b/distributions/normed_gauss1d.hpp
@@ -0,0 +1,58 @@
+#ifndef NGAUSS_MODEL_H_
+#define NGAUSS_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <cmath>
+#include <misc/optvec.hpp>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class normed_gauss1d
+    :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+  {
+  private:
+    normed_gauss1d* do_clone()const
+    {
+      return new normed_gauss1d<T>(*this);
+    }
+
+    const char* do_get_type_name()const
+    {
+      return "1d normed gaussian";
+    }
+  public:
+    normed_gauss1d()
+    {
+      this->push_param_info(param_info<optvec<T> >("x0",0));
+      this->push_param_info(param_info<optvec<T> >("sigma",1));
+    }
+
+    
+  public:
+    optvec<T> do_eval(const optvec<T>& x,const optvec<T>& param)
+    {
+      const double pi=3.14159265358979323846;
+      T x0=get_element(param,0);
+      T sigma=get_element(param,1);
+      if(sigma*sigma<std::numeric_limits<double>::epsilon())
+	{
+	  sigma=std::numeric_limits<double>::epsilon();
+	}
+      T N=1/sqrt(sigma*sigma*pi*2);
+      optvec<T> y=(x-x0)/sigma;
+      return N*exp(-y*y/2.);
+    }
+
+  private:
+    std::string do_get_information()const
+    {
+      return "";
+    }
+  };
+}
+
+
+
+#endif
+//EOF