aboutsummaryrefslogtreecommitdiffstats
path: root/models
diff options
context:
space:
mode:
Diffstat (limited to 'models')
-rw-r--r--models/add_model.hpp169
-rw-r--r--models/beta1d.hpp49
-rw-r--r--models/beta2d.hpp60
-rw-r--r--models/beta2d2.hpp73
-rw-r--r--models/bl1d.hpp56
-rw-r--r--models/bpl1d.hpp56
-rw-r--r--models/bremss.hpp44
-rw-r--r--models/constant.hpp42
-rw-r--r--models/dbeta1d.hpp62
-rw-r--r--models/dbeta2d.hpp92
-rw-r--r--models/dbeta2d2.hpp100
-rw-r--r--models/dbeta2d3.hpp91
-rw-r--r--models/dl_model.hpp178
-rw-r--r--models/dlmodel_template.c36
-rw-r--r--models/func_model.hpp57
-rw-r--r--models/gauss1d.hpp47
-rw-r--r--models/lin1d.hpp42
-rw-r--r--models/models.cc176
-rw-r--r--models/models.hpp36
-rw-r--r--models/mul_model.hpp171
-rw-r--r--models/nbeta1d.hpp49
-rw-r--r--models/nfw1d.hpp45
-rw-r--r--models/pl1d.hpp43
-rw-r--r--models/poly1d.hpp67
-rw-r--r--models/pow_model.hpp120
-rw-r--r--models/strmodel1d.cc82
-rw-r--r--models/strmodel1d.hpp39
-rw-r--r--models/vecn.hpp62
28 files changed, 2144 insertions, 0 deletions
diff --git a/models/add_model.hpp b/models/add_model.hpp
new file mode 100644
index 0000000..ffac9dc
--- /dev/null
+++ b/models/add_model.hpp
@@ -0,0 +1,169 @@
+#ifndef ADD_MODEL_H_
+#define ADD_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename Ty,typename Tx,typename Tp,typename Tstr>
+ class add_model
+ :public model<Ty,Tx,Tp,Tstr>
+ {
+ private:
+ model<Ty,Tx,Tp,Tstr>* do_clone()const
+ {
+ return new add_model<Ty,Tx,Tp,Tstr>(*this);
+ }
+ private:
+ add_model()
+ {
+ }
+
+ private:
+ model<Ty,Tx,Tp,Tstr>* pm1;
+ model<Ty,Tx,Tp,Tstr>* pm2;
+
+ public:
+ add_model(const model<Ty,Tx,Tp,Tstr>& m1,
+ const model<Ty,Tx,Tp,Tstr>& m2)
+ :pm1(m1.clone()),pm2(m2.clone())
+ {
+ int np1=m1.get_num_params();
+ int np2=m2.get_num_params();
+ for(int i=0;i<np1;++i)
+ {
+ param_info<Tp,Tstr> p(m1.get_param_info(i));
+ param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+ this->push_param_info(p1);
+ }
+ for(int i=0;i<np2;++i)
+ {
+ param_info<Tp,Tstr> p(m2.get_param_info(i));
+ param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+ this->push_param_info(p2);
+ }
+ }
+
+ add_model(const add_model& rhs)
+ :pm1(NULL),pm2(NULL)
+ {
+ int np1(0),np2(0);
+ if(rhs.pm1)
+ {
+ pm1=rhs.pm1->clone();
+ np1=rhs.pm1->get_num_params();
+ for(int i=0;i<np1;++i)
+ {
+ param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+ param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+ this->push_param_info(p1);
+ }
+ }
+ if(rhs.pm2)
+ {
+ pm2=rhs.pm2->clone();
+ np2=rhs.pm2->get_num_params();
+ for(int i=0;i<np2;++i)
+ {
+ param_info<Tp,Tstr> p(rhs.pm2->get_param_info(i));
+ param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+ this->push_param_info(p2);
+ }
+ }
+ }
+
+ add_model& operator=(const add_model& rhs)
+ {
+ if(this==&rhs)
+ {
+ return *this;
+ }
+ if(!pm1)
+ {
+ //delete pm1;
+ pm1->destroy();
+ }
+ if(!pm2)
+ {
+ //delete pm2;
+ pm2->destroy();
+ }
+ int np1(0),np2(0);
+ if(rhs.pm1)
+ {
+ pm1=rhs.pm1->clone();
+ np1=rhs.pm1->get_num_params();
+ for(int i=0;i<np1;++i)
+ {
+ param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+ param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+ this->push_param_info(p1);
+ }
+ }
+ if(rhs.pm2)
+ {
+ pm2=rhs.pm2->clone();
+ np2=rhs.pm2->get_num_params();
+ for(int i=0;i<np2;++i)
+ {
+ param_info<Tp,Tstr> p(rhs.pm2->get_param_info(i));
+ param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+ this->push_param_info(p2);
+ }
+ }
+ return *this;
+ }
+
+ ~add_model()
+ {
+ if(!pm1)
+ {
+ //delete pm1;
+ pm1->destroy();
+ }
+ if(!pm2)
+ {
+ //delete pm2;
+ pm2->destroy();
+ }
+ }
+
+ public:
+ Ty do_eval(const Tx& x,const Tp& param)
+ {
+ if(!pm1)
+ {
+ throw opt_exception("incomplete model!");
+ }
+ if(!pm2)
+ {
+ throw opt_exception("incomplete model!");
+ }
+ Tp p1(pm1->get_num_params());
+ Tp p2(pm2->get_num_params());
+ int i=0;
+ int j=0;
+ for(i=0;i<pm1->get_num_params();++i,++j)
+ {
+ set_element(p1,i,get_element(param,j));
+ }
+ for(i=0;i<pm2->get_num_params();++i,++j)
+ {
+ set_element(p2,i,get_element(param,j));
+ }
+ return pm1->eval(x,p1)+pm2->eval(x,p2);
+ }
+ };
+
+ template <typename Ty,typename Tx,typename Tp,typename Tstr>
+ add_model<Ty,Tx,Tp,Tstr> operator+(const model<Ty,Tx,Tp,Tstr>& m1,
+ const model<Ty,Tx,Tp,Tstr>& m2)
+ {
+ return add_model<Ty,Tx,Tp,Tstr>(m1,m2);
+ }
+};
+
+
+
+#endif
+//EOF
diff --git a/models/beta1d.hpp b/models/beta1d.hpp
new file mode 100644
index 0000000..e18aa34
--- /dev/null
+++ b/models/beta1d.hpp
@@ -0,0 +1,49 @@
+#ifndef BETA_MODEL_H_
+#define BETA_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class beta1d
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new beta1d<T>(*this);
+ }
+ public:
+ beta1d()
+ {
+ this->push_param_info(param_info<std::vector<T> >("S0",1));
+ this->push_param_info(param_info<std::vector<T> >("rc",10));
+ this->push_param_info(param_info<std::vector<T> >("beta",2./3.));
+ this->push_param_info(param_info<std::vector<T> >("bkg",0));
+ }
+
+
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ 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 bkg+S0*pow(1+(x*x)/(r_c*r_c),-3*beta+static_cast<T>(.5));
+ }
+
+ std::string do_to_string()const
+ {
+ return "Beta model\n"
+ "S=S0*(1+(r/rc)^2)^(-3*beta+0.5)\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/beta2d.hpp b/models/beta2d.hpp
new file mode 100644
index 0000000..abb6da1
--- /dev/null
+++ b/models/beta2d.hpp
@@ -0,0 +1,60 @@
+#ifndef BETA_MODEL2d_H_
+#define BETA_MODEL2d_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <cassert>
+#include "vecn.hpp"
+
+
+namespace opt_utilities
+{
+
+ template <typename T>
+ class beta2d
+ :public model<T,vecn<T,2>,std::vector<T>,std::string>
+ {
+ private:
+ model<T,vecn<T,2>,std::vector<T> >* do_clone()const
+ {
+ return new beta2d<T>(*this);
+ }
+ public:
+ beta2d()
+ {
+ this->push_param_info(param_info<std::vector<T> >("S0",1));
+ this->push_param_info(param_info<std::vector<T> >("rc1",100));
+ this->push_param_info(param_info<std::vector<T> >("rc2",100));
+ this->push_param_info(param_info<std::vector<T> >("rho",0));
+ this->push_param_info(param_info<std::vector<T> >("x0",100));
+ this->push_param_info(param_info<std::vector<T> >("y0",100));
+ this->push_param_info(param_info<std::vector<T> >("beta",2./3.));
+ this->push_param_info(param_info<std::vector<T> >("bkg",0));
+ }
+
+
+ T do_eval(const vecn<T,2>& xy,const std::vector<T>& param)
+ {
+ T S0=get_element(param,0);
+ T rc1=get_element(param,1);
+ T rc2=get_element(param,2);
+ T rho=get_element(param,3);
+ T x0=get_element(param,4);
+ T y0=get_element(param,5);
+ T beta=get_element(param,6);
+ T bkg=get_element(param,7);
+ T x=xy[0];
+ T y=xy[1];
+
+ T r=(x-x0)*(x-x0)/(rc1*rc1)+(y-y0)*(y-y0)/(rc2*rc2)
+ -2*rho*(x-x0)*(y-y0)/(rc1*rc2);
+ r=r<0?0:r;
+
+ return bkg+S0*pow(1+r,-3*beta+static_cast<T>(.5));
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/beta2d2.hpp b/models/beta2d2.hpp
new file mode 100644
index 0000000..bb473a3
--- /dev/null
+++ b/models/beta2d2.hpp
@@ -0,0 +1,73 @@
+#ifndef BETA_MODEL2d2_H_
+#define BETA_MODEL2d2_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <cassert>
+#include "vecn.hpp"
+
+
+namespace opt_utilities
+{
+
+ template <typename T>
+ class beta2d2
+ :public model<T,vecn<T,2>,std::vector<T>,std::string>
+ {
+ private:
+ model<T,vecn<T,2>,std::vector<T> >* do_clone()const
+ {
+ return new beta2d2<T>(*this);
+ }
+ public:
+ beta2d2()
+ {
+ this->push_param_info(param_info<std::vector<T> >("r0",20));
+ this->push_param_info(param_info<std::vector<T> >("x0",100));
+ this->push_param_info(param_info<std::vector<T> >("y0",100));
+ this->push_param_info(param_info<std::vector<T> >("epsilon",0));
+ this->push_param_info(param_info<std::vector<T> >("theta",100));
+ this->push_param_info(param_info<std::vector<T> >("ampl",3));
+ this->push_param_info(param_info<std::vector<T> >("beta",2./3.));
+ this->push_param_info(param_info<std::vector<T> >("bkg",0));
+ }
+
+
+ T do_eval(const vecn<T,2>& xy,const std::vector<T>& param)
+ {
+ T x=xy[0];
+ T y=xy[1];
+
+ T r0=get_element(param,0);
+ T x0=get_element(param,1);
+ T y0=get_element(param,2);
+ T epsilon=get_element(param,3);
+ T theta=get_element(param,4);
+ T ampl=get_element(param,5);
+ T beta=get_element(param,6);
+ T bkg=get_element(param,7);
+
+ T x_new=(x-x0)*cos(theta)+(y-y0)*sin(theta);
+ T y_new=(y-y0)*cos(theta)-(x-x0)*sin(theta);
+
+ //T _epsilon=sin(epsilon)-0.00001;
+
+ T _epsilon=epsilon;
+
+ // T r=sqrt(x_new*x_new*(1-_epsilon)*(1-_epsilon)
+ // + y_new*y_new);
+
+ //r/=(1-_epsilon);
+ T r_r=x_new*x_new/exp(_epsilon)+y_new*y_new/exp(-_epsilon);
+
+
+ return bkg+ampl*pow(1+r_r/r0/r0,-3*beta+static_cast<T>(.5));
+ }
+
+
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/bl1d.hpp b/models/bl1d.hpp
new file mode 100644
index 0000000..bd0ce7d
--- /dev/null
+++ b/models/bl1d.hpp
@@ -0,0 +1,56 @@
+#ifndef BROKEN_LINE_MODEL_H_
+#define BROKEN_LINE_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class bl1d
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new bl1d<T>(*this);
+ }
+ public:
+ bl1d()
+ {
+ this->push_param_info(param_info<std::vector<T> >("break point y value",1));
+ this->push_param_info(param_info<std::vector<T> >("break point x value",1));
+ this->push_param_info(param_info<std::vector<T> >("slop 1",1));
+ this->push_param_info(param_info<std::vector<T> >("slop 2",1));
+ }
+
+ public:
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ T x_b=get_element(param,0);
+ T f_b=get_element(param,1);
+ T k1=get_element(param,2);
+ T k2=get_element(param,3);
+ if(x<x_b)
+ {
+ return k1*(x-x_b)+f_b;
+ }
+ else
+ {
+ return k2*(x-x_b)+f_b;
+ }
+ }
+
+ private:
+ std::string do_to_string()const
+ {
+ return "broken linear model\n"
+ "y=k1*(x-x_b)+y_b for x<x_b\n"
+ "y=k2*(x-x_b)+y_b otherwise\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/bpl1d.hpp b/models/bpl1d.hpp
new file mode 100644
index 0000000..faf1336
--- /dev/null
+++ b/models/bpl1d.hpp
@@ -0,0 +1,56 @@
+#ifndef BROKEN_POWER_LAW_MODEL_H_
+#define BROKEN_POWER_LAW_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class bpl1d
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new bpl1d<T>(*this);
+ }
+ public:
+ bpl1d()
+ {
+ this->push_param_info(param_info<std::vector<T> >("bpx",1));
+ this->push_param_info(param_info<std::vector<T> >("bpy",1));
+ this->push_param_info(param_info<std::vector<T> >("gamma1",1));
+ this->push_param_info(param_info<std::vector<T> >("gamma2",1));
+ }
+
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ T x_b=get_element(param,0);
+ T f_b=get_element(param,1);
+ T gamma1=get_element(param,2);
+ T gamma2=get_element(param,3);
+ if(x<x_b)
+ {
+ return f_b*pow(x,gamma1)/pow(x_b,gamma1);
+ }
+ else
+ {
+ return f_b*pow(x,gamma2)/pow(x_b,gamma2);
+ }
+ }
+
+
+ private:
+ std::string do_to_string()const
+ {
+ return "broken power law\n"
+ "y=y_b*(x/x_b)^gamma1 for x<x_b\n"
+ "y=y_b*(x/x_b)^gamma2 otherwise\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/bremss.hpp b/models/bremss.hpp
new file mode 100644
index 0000000..9956a3e
--- /dev/null
+++ b/models/bremss.hpp
@@ -0,0 +1,44 @@
+#ifndef BREMSS_MODEL_H_
+#define BREMSS_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class bremss
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new bremss<T>(*this);
+ }
+ public:
+ bremss()
+ {
+ this->push_param_info(param_info<std::vector<T> >("norm",1));
+ this->push_param_info(param_info<std::vector<T> >("kT",1));
+ }
+
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ T norm=get_element(param,0);
+ T kT=get_element(param,1);
+
+ return norm*sqrt(kT)*exp(-x/kT);
+ }
+
+ private:
+ std::string do_to_string()const
+ {
+ return "Simplified bremss model\n"
+ "flux=norm*kT^0.5*e^{-E/kT}\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/constant.hpp b/models/constant.hpp
new file mode 100644
index 0000000..8d06279
--- /dev/null
+++ b/models/constant.hpp
@@ -0,0 +1,42 @@
+#ifndef CONSTANT_MODEL_H_
+#define CONSTANT_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class constant
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new constant<T>(*this);
+ }
+ public:
+ constant()
+ {
+ this->push_param_info(param_info<std::vector<T> >("c",1));
+ }
+
+ public:
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ //return x*param[0]+param[1];
+ return get_element(param,0);
+ }
+
+ private:
+ std::string do_to_string()const
+ {
+ return "Constant\n"
+ "y=C\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/dbeta1d.hpp b/models/dbeta1d.hpp
new file mode 100644
index 0000000..71a1f52
--- /dev/null
+++ b/models/dbeta1d.hpp
@@ -0,0 +1,62 @@
+#ifndef DBETA_MODEL_H_
+#define DBETA_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+namespace opt_utilities
+{
+ template <typename T>
+ class dbeta1d
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new dbeta1d<T>(*this);
+ }
+ public:
+ dbeta1d()
+ {
+ this->push_param_info(param_info<std::vector<T> >("S01",1));
+ this->push_param_info(param_info<std::vector<T> >("rc1",30));
+ this->push_param_info(param_info<std::vector<T> >("beta1",.6));
+
+ this->push_param_info(param_info<std::vector<T> >("S02",.5));
+ this->push_param_info(param_info<std::vector<T> >("rc2",20));
+ this->push_param_info(param_info<std::vector<T> >("beta2",.4));
+
+
+ this->push_param_info(param_info<std::vector<T> >("bkg",0));
+ }
+
+
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ T S01=get_element(param,0);
+ T r_c1=get_element(param,1);
+ T beta1=get_element(param,2);
+
+ T S02=get_element(param,3);
+ T r_c2=get_element(param,4);
+ T beta2=get_element(param,5);
+
+
+ T bkg=get_element(param,6);
+
+ return bkg+S01*pow(1+(x*x)/(r_c1*r_c1),-3*beta1+static_cast<T>(.5))
+ +S02*pow(1+(x*x)/(r_c2*r_c2),-3*beta2+static_cast<T>(.5));
+ }
+
+ private:
+ std::string do_to_string()const
+ {
+ return "double 1d beta model\n"
+ "S=beta(S01,beta1,rc1)+beta(S02,beta2,rc2)\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/dbeta2d.hpp b/models/dbeta2d.hpp
new file mode 100644
index 0000000..6811270
--- /dev/null
+++ b/models/dbeta2d.hpp
@@ -0,0 +1,92 @@
+#ifndef DBETA_MODEL2d_H_
+#define DBETA_MODEL2d_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <cassert>
+#include "vecn.hpp"
+
+namespace opt_utilities
+{
+
+ template <typename T>
+ class dbeta2d
+ :public model<T,vecn<T,2>,std::vector<T>,std::string>
+ {
+ private:
+ model<T,vecn<T,2>,std::vector<T> >* do_clone()const
+ {
+ return new dbeta2d<T>(*this);
+ }
+
+
+
+ public:
+ dbeta2d()
+ {
+ push_param_info(param_info<std::vector<T> >("S01",1));
+ push_param_info(param_info<std::vector<T> >("rc11",50));
+ push_param_info(param_info<std::vector<T> >("rc21",50));
+ push_param_info(param_info<std::vector<T> >("rho1",0));
+ push_param_info(param_info<std::vector<T> >("x01",200));
+ push_param_info(param_info<std::vector<T> >("y01",200));
+ push_param_info(param_info<std::vector<T> >("beta1",2./3.));
+
+ push_param_info(param_info<std::vector<T> >("S02",1));
+ push_param_info(param_info<std::vector<T> >("rc12",60));
+ push_param_info(param_info<std::vector<T> >("rc22",60));
+ push_param_info(param_info<std::vector<T> >("rho2",0));
+ push_param_info(param_info<std::vector<T> >("x02",200));
+ push_param_info(param_info<std::vector<T> >("y02",200));
+ push_param_info(param_info<std::vector<T> >("beta2",2./2.5));
+ push_param_info(param_info<std::vector<T> >("bkg",0));
+ }
+
+ T do_eval(const vecn<T,2>& xy,const std::vector<T>& param)
+ {
+ T S01=get_element(param,0);
+ T rc11=get_element(param,1);
+ T rc21=get_element(param,2);
+ T rho1=get_element(param,3);
+ T x01=get_element(param,4);
+ T y01=get_element(param,5);
+ T beta1=get_element(param,6);
+ T S02=get_element(param,7);
+ T rc12=get_element(param,8);
+ T rc22=get_element(param,9);
+ T rho2=get_element(param,10);
+ T x02=get_element(param,11);
+ T y02=get_element(param,12);
+ T beta2=get_element(param,13);
+ T bkg=get_element(param,14);
+
+
+ T x=xy[0];
+ T y=xy[1];
+
+ rho1=rho1>1?1:rho1;
+ rho1=rho1<-1?-1:rho1;
+ rho2=rho2>1?1:rho2;
+ rho2=rho2<-1?-1:rho2;
+
+ T r1=(x-x01)*(x-x01)/(rc11*rc11)+(y-y01)*(y-y01)/(rc21*rc21)
+ -2*rho1*(x-x01)*(y-y01)/(rc11*rc21);
+
+ T r2=(x-x02)*(x-x02)/(rc12*rc12)+(y-y02)*(y-y02)/(rc22*rc22)
+ -2*rho2*(x-x02)*(y-y02)/(rc12*rc22);
+ // r1=r1<0?0:r1;
+ //r2=r2<0?0:r2;
+ assert(r1>=0);
+ assert(r2>=0);
+
+
+ return bkg+S01*pow(1+r1,-3*beta1+static_cast<T>(.5))+
+ S02*pow(1+r2,-3*beta2+static_cast<T>(.5));
+ }
+ };
+};
+
+
+
+#endif
+//EOF
+
diff --git a/models/dbeta2d2.hpp b/models/dbeta2d2.hpp
new file mode 100644
index 0000000..d968dbc
--- /dev/null
+++ b/models/dbeta2d2.hpp
@@ -0,0 +1,100 @@
+#ifndef DDBETA_OOOMODEL2d2_H_
+#define DDBETA_OOOMODEL2d2_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <cassert>
+#include "vecn.hpp"
+
+
+namespace opt_utilities
+{
+
+ template <typename T>
+ class dbeta2d2
+ :public model<T,vecn<T,2>,std::vector<T>,std::string>
+ {
+ private:
+ model<T,vecn<T,2>,std::vector<T> >* do_clone()const
+ {
+ return new dbeta2d2<T>(*this);
+ }
+ public:
+ dbeta2d2()
+ {
+ this->push_param_info(param_info<std::vector<T> >("ampl1",1));//0
+ this->push_param_info(param_info<std::vector<T> >("r01",1));//1
+ this->push_param_info(param_info<std::vector<T> >("x01",100));//2
+ this->push_param_info(param_info<std::vector<T> >("y01",100));//3
+ this->push_param_info(param_info<std::vector<T> >("theta1",0));//4
+ this->push_param_info(param_info<std::vector<T> >("beta1",2./3.));//5
+ this->push_param_info(param_info<std::vector<T> >("epsilon1",0));//6
+
+
+ this->push_param_info(param_info<std::vector<T> >("ampl2",1));//7
+ this->push_param_info(param_info<std::vector<T> >("r02",1));//8
+ this->push_param_info(param_info<std::vector<T> >("x02",100));//9
+ this->push_param_info(param_info<std::vector<T> >("y02",100));//10
+ this->push_param_info(param_info<std::vector<T> >("theta2",0));//11
+ this->push_param_info(param_info<std::vector<T> >("beta2",2./4.));//12
+ this->push_param_info(param_info<std::vector<T> >("epsilon2",0));//13
+
+ this->push_param_info(param_info<std::vector<T> >("bkg",0));//14
+ }
+
+
+ T do_eval(const vecn<T,2>& xy,const std::vector<T>& param)
+ {
+ T x=xy[0];
+ T y=xy[1];
+
+ T ampl1=get_element(param,0);
+ T r01=get_element(param,1);
+ T x01=get_element(param,2);
+ T y01=get_element(param,3);
+
+ T theta1=get_element(param,4);
+
+ T beta1=get_element(param,5);
+ T epsilon1=get_element(param,6);
+
+ T ampl2=get_element(param,7);
+ T r02=get_element(param,8);
+ T x02=get_element(param,9);
+ T y02=get_element(param,10);
+
+ T theta2=get_element(param,11);
+
+ T beta2=get_element(param,12);
+
+ T epsilon2=get_element(param,13);
+ T bkg=get_element(param,14);
+
+ T x_new1=(x-x01)*cos(theta1)+(y-y01)*sin(theta1);
+ T y_new1=(y-y01)*cos(theta1)-(x-x01)*sin(theta1);
+
+ //T r1=sqrt(x_new1*x_new1*(1-epsilon1)*(1-epsilon1)
+ // + y_new1*y_new1);
+ //r1/=(1-epsilon1);
+
+ T r1_r1=x_new1*x_new1/exp(epsilon1/30)+y_new1*y_new1/exp(-epsilon1/30);
+
+ T x_new2=(x-x02)*cos(theta2)+(y-y02)*sin(theta2);
+ T y_new2=(y-y02)*cos(theta2)-(x-x02)*sin(theta2);
+
+ T r2_r2=x_new2*x_new2/exp(epsilon2/30)+y_new2*y_new2/exp(-epsilon2/30);
+ // T r2=sqrt(x_new2*x_new2*(1-epsilon2)*(1-epsilon2)
+ // + y_new2*y_new2);
+ //r2/=(1-epsilon2);
+
+
+
+ return bkg+ampl1*pow(1+(r1_r1/r01/r01),-3*beta1+static_cast<T>(.5))
+ +ampl2*pow(1+(r2_r2/r02/r02),-3*beta2+static_cast<T>(.5));
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/dbeta2d3.hpp b/models/dbeta2d3.hpp
new file mode 100644
index 0000000..e2edb78
--- /dev/null
+++ b/models/dbeta2d3.hpp
@@ -0,0 +1,91 @@
+#ifndef DDBETA3_MODEL2d2_H_
+#define DDBETA3_MODEL2d2_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <cassert>
+#include "vecn.hpp"
+
+
+namespace opt_utilities
+{
+
+ template <typename T>
+ class dbeta2d3
+ :public model<T,vecn<T,2>,std::vector<T>,std::string>
+ {
+ private:
+ model<T,vecn<T,2>,std::vector<T> >* do_clone()const
+ {
+ return new dbeta2d3<T>(*this);
+ }
+ public:
+ dbeta2d3()
+ {
+
+ this->push_param_info(param_info<std::vector<T> >("x0",256));//1
+ this->push_param_info(param_info<std::vector<T> >("y0",256));//2
+ this->push_param_info(param_info<std::vector<T> >("epsilon",0));//3
+ this->push_param_info(param_info<std::vector<T> >("theta",0));//4
+
+ this->push_param_info(param_info<std::vector<T> >("ampl1",1));//5
+ this->push_param_info(param_info<std::vector<T> >("beta1",0.6));//6
+ this->push_param_info(param_info<std::vector<T> >("r01",30));//7
+
+ this->push_param_info(param_info<std::vector<T> >("ampl2",.5));//8
+ this->push_param_info(param_info<std::vector<T> >("beta2",.4));//9
+ this->push_param_info(param_info<std::vector<T> >("r02",20));//10
+
+ this->push_param_info(param_info<std::vector<T> >("bkg",0));//11
+
+
+ }
+
+
+ T do_eval(const vecn<T,2>& xy,const std::vector<T>& param)
+ {
+ T x=xy[0];
+ T y=xy[1];
+
+ T x0=get_element(param,0);
+ T y0=get_element(param,1);
+ T epsilon=get_element(param,2);
+ T theta=get_element(param,3);
+
+ T ampl1=(get_element(param,4));
+ T beta1=(get_element(param,5));
+ T r01=(get_element(param,6));
+
+ T ampl2=(get_element(param,7));
+ T beta2=(get_element(param,8));
+ T r02=(get_element(param,9));
+
+ T bkg=get_element(param,10);
+
+ T x_new1=(x-x0)*cos(theta)+(y-y0)*sin(theta);
+ T y_new1=(y-y0)*cos(theta)-(x-x0)*sin(theta);
+
+
+ T r1_r1=x_new1*x_new1/exp(epsilon/30)+y_new1*y_new1/exp(-epsilon/30);
+ //T r1=sqrt(x_new1*x_new1*(1-epsilon)*(1-epsilon)+y_new1*y_new1)/(1-epsilon);
+
+ T x_new2=(x-x0)*cos(theta)+(y-y0)*sin(theta);
+ T y_new2=(y-y0)*cos(theta)-(x-x0)*sin(theta);
+
+ T r2_r2=x_new2*x_new2/exp(epsilon/30)+y_new2*y_new2/exp(-epsilon/30);
+ //T r2=sqrt(x_new2*x_new2*(1-epsilon)*(1-epsilon)+y_new2*y_new2)/(1-epsilon);
+
+
+
+ return bkg+ampl1*pow(1+(r1_r1/r01/r01),-3*beta1+static_cast<T>(.5))
+ +ampl2*pow(1+(r2_r2/r02/r02),-3*beta2+static_cast<T>(.5));
+
+ //return bkg+pow(1+r1*r1/r01/r01,-3*beta1+static_cast<T>(.5))+
+ //pow(1+r2*r2/r02/r02,-3*beta2+static_cast<T>(.5));
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/dl_model.hpp b/models/dl_model.hpp
new file mode 100644
index 0000000..2fcf227
--- /dev/null
+++ b/models/dl_model.hpp
@@ -0,0 +1,178 @@
+#ifdef __linux__
+
+#ifndef DL_MODEL_H_
+#define DL_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+#include <string>
+#include <sstream>
+#include <dlfcn.h>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class dl_model
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ T (*calc_model)(T x,const T* p);
+ int nparams;
+ mutable void* handle;
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ dl_model<T>* result=new dl_model<T>(*this);
+ this->handle=NULL;
+ return result;
+ }
+
+ // public:
+ public:
+ dl_model()
+ :handle(NULL)
+ {}
+
+
+ public:
+ dl_model(const char* file_name)
+ :handle(NULL)
+ {
+
+ handle=dlopen(file_name,RTLD_LAZY);
+
+ if(!handle)
+ {
+ throw opt_exception("faild loading object");
+ }
+
+ calc_model=(T (*)(T,const T*))dlsym(handle,"calc_model");
+
+ if(!calc_model)
+ {
+ throw opt_exception("symble undefined");
+ }
+
+ const char* (*get_param_name)(int)
+ =(const char* (*)(int))dlsym(handle,"get_param_name");
+
+ if(!get_param_name)
+ {
+ throw opt_exception("symble undefined");
+ }
+
+ int (*get_num_params)()
+ =(int (*)())dlsym(handle,"get_num_params");
+
+ if(!get_num_params)
+ {
+ throw opt_exception("symble undefined");
+ if(!get_num_params)
+ {
+ throw opt_exception("symble undefined");
+ } }
+
+ T (*get_default_value)(int)
+ =(T (*)(int))dlsym(handle,"get_default_value");
+
+ if(!get_default_value)
+ {
+ throw opt_exception("symble undefined");
+ }
+
+ nparams=get_num_params();
+
+ for(int i=0;i!=nparams;++i)
+ {
+ this->push_param_info(param_info<std::vector<T> >(get_param_name(i),
+ get_default_value(i)));
+ }
+ }
+
+ ~dl_model()
+ {
+ if(handle)
+ {
+ dlclose(handle);
+ }
+ }
+
+ void bind(const char* file_name)
+ {
+ if(handle)
+ {
+ dlclose(handle);
+ }
+ this->clear_param_info();
+ handle=dlopen(file_name,RTLD_LAZY);
+
+ if(!handle)
+ {
+ throw opt_exception("faild loading object");
+ }
+
+ calc_model=(T (*)(T,const T*))dlsym(handle,"calc_model");
+
+ if(!calc_model)
+ {
+ throw opt_exception("symble undefined");
+ }
+
+ const char* (*get_param_name)(int)
+ =(const char* (*)(int))dlsym(handle,"get_param_name");
+
+ if(!get_param_name)
+ {
+ throw opt_exception("symble undefined");
+ }
+
+
+ int (*get_num_params)()
+ =(int (*)())dlsym(handle,"get_num_params");
+
+ if(!get_num_params)
+ {
+ throw opt_exception("symble undefined");
+ }
+
+
+ T (*get_default_value)(int)
+ =(T (*)(int))dlsym(handle,"get_default_value");
+
+
+ if(!get_default_value)
+ {
+ throw opt_exception("symble undefined");
+ }
+
+
+ nparams=get_num_params();
+ for(int i=0;i!=nparams;++i)
+ {
+ this->push_param_info(param_info<std::vector<T> >(get_param_name(i),
+ get_default_value(i)));
+ }
+ }
+
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ if(handle==NULL)
+ {
+ throw opt_exception("dl object unloaded");
+ }
+ return calc_model(x,&get_element(param,0));
+ }
+
+ std::string do_to_string()const
+ {
+ return "Dynamical load model\n"
+ "Should be loaded from an shared object file\n";
+ }
+ };
+};
+
+
+
+#endif
+#endif
+//EOF
diff --git a/models/dlmodel_template.c b/models/dlmodel_template.c
new file mode 100644
index 0000000..fd7aadb
--- /dev/null
+++ b/models/dlmodel_template.c
@@ -0,0 +1,36 @@
+#include <math.h>
+
+char p1name[2]="k";
+char p2name[2]="b";
+
+int main(int argc,char* argv[])
+{}
+
+
+int get_num_params()
+{
+ return 2;
+}
+
+const char* get_param_name(int n)
+{
+ if(n==0)
+ {
+ return p1name;
+ }
+ return p2name;
+}
+
+double get_default_value(int n)
+{
+ if(n==0)
+ return 1;
+ return 0;
+}
+
+double calc_model(double x,double p[])
+{
+ return p[0]*x+p[1];
+}
+
+
diff --git a/models/func_model.hpp b/models/func_model.hpp
new file mode 100644
index 0000000..4eda8d5
--- /dev/null
+++ b/models/func_model.hpp
@@ -0,0 +1,57 @@
+#ifndef FUNC_MODEL_H_
+#define FUNC_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+#include <string>
+#include <sstream>
+namespace opt_utilities
+{
+ template <typename T>
+ class func_model
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ T (*func)(T x,const T* const& p);
+ int nparams;
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new func_model<T>(*this);
+ }
+
+ // public:
+ private:
+ func_model()
+ {}
+
+
+ public:
+ func_model(T (*_func)(T x,const T* const& p),int n)
+ :func(_func),nparams(n)
+ {
+ for(int i=0;i!=n;++i)
+ {
+ std::ostringstream oss;
+ oss<<i;
+ this->push_param_info(param_info<std::vector<T> >(oss.str(),0));
+ }
+ }
+
+
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ return func(x,&get_element(param,0));
+ }
+ private:
+ std::string do_to_string()const
+ {
+ return "Wrapper for necked C function\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/gauss1d.hpp b/models/gauss1d.hpp
new file mode 100644
index 0000000..958cec9
--- /dev/null
+++ b/models/gauss1d.hpp
@@ -0,0 +1,47 @@
+#ifndef GAUSS_MODEL_H_
+#define GAUSS_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class gauss1d
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new gauss1d<T>(*this);
+ }
+ public:
+ gauss1d()
+ {
+ this->push_param_info(param_info<std::vector<T> >("N",1));
+ this->push_param_info(param_info<std::vector<T> >("x0",0));
+ this->push_param_info(param_info<std::vector<T> >("sigma",1));
+ }
+
+ public:
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ T N=get_element(param,0);
+ T x0=get_element(param,1);
+ T sigma=get_element(param,2);
+ T y=(x-x0)/sigma;
+ return N*exp(-y*y);
+ }
+
+ private:
+ std::string do_to_string()const
+ {
+ return "Gaussian model\n"
+ "y=N*exp(-(x-x0)^2/sigma^2\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/lin1d.hpp b/models/lin1d.hpp
new file mode 100644
index 0000000..b0ddd31
--- /dev/null
+++ b/models/lin1d.hpp
@@ -0,0 +1,42 @@
+#ifndef LINEAR_MODEL_H_
+#define LINEAR_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class lin1d
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new lin1d<T>(*this);
+ }
+ public:
+ lin1d()
+ {
+ this->push_param_info(param_info<std::vector<T> >("k",1));
+ this->push_param_info(param_info<std::vector<T> >("b",0));
+ }
+
+ public:
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ return x*get_element(param,0)+get_element(param,1);
+ }
+
+ private:
+ std::string do_to_string()const
+ {
+ return "linear model\n"
+ "y=k*x+b\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/models.cc b/models/models.cc
new file mode 100644
index 0000000..7481b88
--- /dev/null
+++ b/models/models.cc
@@ -0,0 +1,176 @@
+#include "models.hpp"
+#include <core/opt_exception.hpp>
+#include "gauss1d.hpp"
+#include "bl1d.hpp"
+#include "nfw1d.hpp"
+#include "bpl1d.hpp"
+#include "beta1d.hpp"
+#include "nbeta1d.hpp"
+#include "dbeta1d.hpp"
+#include "lin1d.hpp"
+#include "pl1d.hpp"
+#include "poly1d.hpp"
+#include "bremss.hpp"
+#include "beta2d2.hpp"
+#include "beta2d.hpp"
+#include "dbeta2d2.hpp"
+#include "dbeta2d3.hpp"
+#include "dbeta2d.hpp"
+#include "dl_model.hpp"
+#include <iostream>
+using namespace std;
+
+
+namespace opt_utilities
+{
+ strmodel1d strm1d;
+ std::map<std::string,model<double,double,std::vector<double>,std::string>* > model_map;
+ std::map<std::string,model<double,vecn<double,2>,std::vector<double>,std::string>* > model2d_map;
+ std::list<std::string> get_model_name_list()
+ {
+ std::list<std::string> result;
+ for(std::map<std::string,model<double,double,std::vector<double>,std::string>* >::iterator i=model_map.begin();
+ i!=model_map.end();++i)
+ {
+ result.push_back(i->first);
+ }
+ return result;
+ }
+
+ std::list<std::string> get_model2d_name_list()
+ {
+ std::list<std::string> result;
+ for(map<std::string,model<double,vecn<double,2>,std::vector<double>,std::string>* > ::iterator i=model2d_map.begin();
+ i!=model2d_map.end();++i)
+ {
+ result.push_back(i->first);
+ }
+ return result;
+ }
+
+model<double,double,std::vector<double>,std::string>& get_1dmodel_by_name(const char* name)
+ {
+ std::map<std::string,model<double,double,std::vector<double>,std::string >* >::iterator iter;
+ iter=model_map.find(name);
+ if(iter==model_map.end()||iter->second==0)
+ {
+ throw opt_exception("model does not exist");
+ }
+ return *(iter->second);
+ }
+
+ strmodel1d& get_strm1d()
+ {
+ return strm1d;
+ }
+
+model<double,vecn<double,2>,std::vector<double>,std::string>& get_2dmodel_by_name(const char* name)
+ {
+ std::map<std::string,model<double,vecn<double,2>,std::vector<double>,std::string>* >::iterator iter;
+ iter=model2d_map.find(name);
+ if(iter==model2d_map.end()||iter->second==0)
+ {
+ throw opt_exception("model does not exist");
+ }
+ return *(iter->second);
+ }
+
+ int get_n_1dmodels()
+ {
+ return model_map.size();
+ }
+
+ int get_n_2dmodels()
+ {
+ return model2d_map.size();
+ }
+
+ class model_map_keeper_class
+ {
+ private:
+ void init_model_map()
+ {
+
+ //#define DECL_POLY(n) model_map["poly1d##n"]=new poly1d<double,n>;
+
+
+ model_map["lin1d"]=new lin1d<double>;
+ model_map["pl1d"]=new pl1d<double>;
+ model_map["bl1d"]=new bl1d<double>;
+ model_map["bpl1d"]=new bpl1d<double>;
+ model_map["beta1d"]=new beta1d<double>;
+ model_map["bremss"]=new bremss<double>;
+ model_map["nbeta1d"]=new nbeta1d<double>;
+ model_map["2beta1d"]=new dbeta1d<double>;
+ model_map["nfw"]=new nfw1d<double>;
+ model_map["gauss1d"]=new gauss1d<double>;
+ model_map["poly1d2"]=new poly1d<double,2>;
+ model_map["poly1d3"]=new poly1d<double,3>;
+ model_map["poly1d4"]=new poly1d<double,4>;
+ model_map["poly1d5"]=new poly1d<double,5>;
+#ifdef __linux__
+ model_map["dlmodel"]=new dl_model<double>;
+#endif
+ //DECL_POLY(7)
+ }
+
+ void release_model_map()
+ {
+ for(std::map<std::string,model<double,double,std::vector<double>,std::string>* >::iterator i=model_map.begin();
+ i!=model_map.end();++i)
+ {
+ delete i->second;
+ }
+ }
+ public:
+ model_map_keeper_class()
+ {
+ init_model_map();
+ std::cerr<<"1d models Initialized"<<std::endl;
+ }
+
+ ~model_map_keeper_class()
+ {
+ release_model_map();
+ std::cerr<<"1d models Released"<<std::endl;
+ }
+
+ }model_map_keeper;
+
+ class model2d_map_keeper_class
+ {
+ private:
+ void init_model_map()
+ {
+ model2d_map["beta2d"]=new beta2d<double>;
+ model2d_map["beta2d2"]=new beta2d2<double>;
+ model2d_map["dbeta2d"]=new dbeta2d<double>;
+ model2d_map["dbeta2d2"]=new dbeta2d2<double>;
+ model2d_map["dbeta2d3"]=new dbeta2d3<double>;
+ }
+
+ void release_model_map()
+ {
+ for(std::map<std::string,model<double,vecn<double,2>,std::vector<double>,std::string>* >::iterator i=model2d_map.begin();
+ i!=model2d_map.end();++i)
+ {
+ delete i->second;
+ }
+ }
+ public:
+ model2d_map_keeper_class()
+ {
+ init_model_map();
+ std::cerr<<"2d models Initialized"<<std::endl;
+ }
+
+ ~model2d_map_keeper_class()
+ {
+ release_model_map();
+ std::cerr<<"2d models Released"<<std::endl;
+ }
+
+ }model2d_map_keeper;
+
+
+};
diff --git a/models/models.hpp b/models/models.hpp
new file mode 100644
index 0000000..79201e8
--- /dev/null
+++ b/models/models.hpp
@@ -0,0 +1,36 @@
+#ifndef MODELS_HPP
+#define MODELS_HPP
+
+#include <core/fitter.hpp>
+#include <map>
+#include <string>
+#include <list>
+#include "vecn.hpp"
+#include "strmodel1d.hpp"
+
+
+
+namespace opt_utilities
+{
+ extern std::map<std::string,model<double,double,std::vector<double>,std::string>* > model_map;
+ extern std::map<std::string,model<double,vecn<double,2>,std::vector<double>,std::string >* > model2d_map;
+
+ extern strmodel1d strm1d;
+ extern std::list<std::string> get_model_name_list();
+ extern int get_n_1dmodels();
+ // extern void init_model_map();
+ // extern void release_model_map();
+
+ extern std::list<std::string> get_model2d_name_list();
+ // extern void init_model2d_map();
+ // extern void release_model2d_map();
+ extern int get_n_2dmodels();
+
+ extern model<double,double,std::vector<double>,std::string >& get_1dmodel_by_name(const char*);
+ extern model<double,vecn<double,2>,std::vector<double>,std::string >& get_2dmodel_by_name(const char*);
+
+ extern strmodel1d& get_strm1d();
+}
+
+
+#endif
diff --git a/models/mul_model.hpp b/models/mul_model.hpp
new file mode 100644
index 0000000..0e64c1b
--- /dev/null
+++ b/models/mul_model.hpp
@@ -0,0 +1,171 @@
+#ifndef MUL_MODEL_H_
+#define MUL_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename Ty,typename Tx,typename Tp,typename Tstr>
+ class mul_model
+ :public model<Ty,Tx,Tp,Tstr>
+ {
+ private:
+ model<Ty,Tx,Tp,Tstr>* do_clone()const
+ {
+ return new mul_model<Ty,Tx,Tp,Tstr>(*this);
+ }
+ private:
+ mul_model()
+ {
+ }
+
+ private:
+ model<Ty,Tx,Tp,Tstr>* pm1;
+ model<Ty,Tx,Tp,Tstr>* pm2;
+
+ public:
+ mul_model(const model<Ty,Tx,Tp,Tstr>& m1,
+ const model<Ty,Tx,Tp,Tstr>& m2)
+ :pm1(m1.clone()),pm2(m2.clone())
+ {
+ int np1=m1.get_num_params();
+ int np2=m2.get_num_params();
+ for(int i=0;i<np1;++i)
+ {
+ param_info<Tp,Tstr> p(m1.get_param_info(i));
+ param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+ this->push_param_info(p1);
+ }
+ for(int i=0;i<np2;++i)
+ {
+ param_info<Tp,Tstr> p(m2.get_param_info(i));
+ param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+ this->push_param_info(p2);
+ }
+ }
+
+ mul_model(const mul_model& rhs)
+ :pm1(NULL),pm2(NULL)
+ {
+ int np1(0),np2(0);
+ if(rhs.pm1)
+ {
+ pm1=rhs.pm1->clone();
+ np1=rhs.pm1->get_num_params();
+ for(int i=0;i<np1;++i)
+ {
+ param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+ param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+ this->push_param_info(p1);
+ }
+ }
+ if(rhs.pm2)
+ {
+ pm2=rhs.pm2->clone();
+ np2=rhs.pm2->get_num_params();
+ for(int i=0;i<np2;++i)
+ {
+ param_info<Tp,Tstr> p(rhs.pm2->get_param_info(i));
+ param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+ this->push_param_info(p2);
+ }
+ }
+ }
+
+ mul_model& operator=(const mul_model& rhs)
+ {
+ if(this==&rhs)
+ {
+ return *this;
+ }
+ if(!pm1)
+ {
+ //delete pm1;
+ pm1->destroy();
+ }
+ if(!pm2)
+ {
+ //delete pm2;
+ pm2->destroy();
+ }
+ int np1(0),np2(0);
+ if(rhs.pm1)
+ {
+ pm1=rhs.pm1->clone();
+ np1=rhs.pm1->get_num_params();
+ for(int i=0;i<np1;++i)
+ {
+ param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+ param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+ this->push_param_info(p1);
+ }
+ }
+ if(rhs.pm2)
+ {
+ pm2=rhs.pm2->clone();
+ np2=rhs.pm2->get_num_params();
+ for(int i=0;i<np2;++i)
+ {
+ param_info<Tp,Tstr> p(rhs.pm2->get_param_info(i));
+ param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+ this->push_param_info(p2);
+ }
+ }
+ return *this;
+ }
+
+
+ ~mul_model()
+ {
+ if(!pm1)
+ {
+ //delete pm1;
+ pm1->destroy();
+ }
+ if(!pm2)
+ {
+ //delete pm2;
+ pm2->destroy();
+ }
+ }
+
+ public:
+ Ty do_eval(const Tx& x,const Tp& param)
+ {
+ if(!pm1)
+ {
+ throw opt_exception("incomplete model!");
+ }
+ if(!pm2)
+ {
+ throw opt_exception("incomplete model!");
+ }
+ Tp p1(pm1->get_num_params());
+ Tp p2(pm2->get_num_params());
+ int i=0;
+ int j=0;
+ for(i=0;i<pm1->get_num_params();++i,++j)
+ {
+ set_element(p1,i,get_element(param,j));
+ }
+ for(i=0;i<pm2->get_num_params();++i,++j)
+ {
+ set_element(p2,i,get_element(param,j));
+ }
+ return pm1->eval(x,p1)*pm2->eval(x,p2);
+ }
+ };
+
+ template <typename Ty,typename Tx,typename Tp>
+ mul_model<Ty,Tx,Tp,Tstr> operator*(const model<Ty,Tx,Tp,Tstr>& m1,
+ const model<Ty,Tx,Tp,Tstr>& m2)
+ {
+ return mul_model<Ty,Tx,Tp,Tstr>(m1,m2);
+ }
+
+};
+
+
+
+#endif
+//EOF
diff --git a/models/nbeta1d.hpp b/models/nbeta1d.hpp
new file mode 100644
index 0000000..3806fe9
--- /dev/null
+++ b/models/nbeta1d.hpp
@@ -0,0 +1,49 @@
+#ifndef NBETA_MODEL_H_
+#define NBETA_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+namespace opt_utilities
+{
+ template <typename T>
+ class nbeta1d
+ :public model<T,T,std::vector<T> ,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new nbeta1d<T>(*this);
+ }
+ public:
+ nbeta1d()
+ {
+ this->push_param_info(param_info<std::vector<T> >("S0",1));
+ this->push_param_info(param_info<std::vector<T> >("rc",10));
+ this->push_param_info(param_info<std::vector<T> >("beta",2./3.));
+ this->push_param_info(param_info<std::vector<T> >("bkg",0));
+ }
+
+
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ 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 bkg+S0*pow(1+(x*x)/(r_c*r_c),-3./2.*beta);
+ }
+
+ private:
+ std::string do_to_string()const
+ {
+ return "density beta model\n"
+ "n=n0*(1+(r/rc)^2)^(-1.5*beta)\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/nfw1d.hpp b/models/nfw1d.hpp
new file mode 100644
index 0000000..9c2dfd4
--- /dev/null
+++ b/models/nfw1d.hpp
@@ -0,0 +1,45 @@
+#ifndef NFW_MODEL_H_
+#define NFW_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+namespace opt_utilities
+{
+ template <typename T>
+ class nfw1d
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new nfw1d<T>(*this);
+ }
+ public:
+ nfw1d()
+ {
+ this->push_param_info(param_info<std::vector<T> >("rho0",1.));
+ this->push_param_info(param_info<std::vector<T> >("rs",1.));
+ }
+
+
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ T rho0=get_element(param,0);
+ T rs=get_element(param,1);
+
+
+ return rho0/(x/rs*(1+x/rs)*(1+x/rs));
+ }
+
+ std::string do_to_string()const
+ {
+ return "NFW model\n"
+ "y=rho0/(r/rs*(1+r/rs)^2\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/pl1d.hpp b/models/pl1d.hpp
new file mode 100644
index 0000000..50b4c50
--- /dev/null
+++ b/models/pl1d.hpp
@@ -0,0 +1,43 @@
+#ifndef POWER_LAW_MODEL_H_
+#define POWER_LAW_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class pl1d
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new pl1d<T>(*this);
+ }
+ public:
+ pl1d()
+ {
+ this->push_param_info(param_info<std::vector<T> >("Ampl",1));
+ this->push_param_info(param_info<std::vector<T> >("gamma",1));
+ }
+
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ T A=get_element(param,0);
+ T gamma=get_element(param,1);
+ return A*pow(x,gamma);
+ }
+
+ private:
+ std::string do_to_string()const
+ {
+ return "Simple power law model\n"
+ "y=A*x^gamma\n";
+ }
+ };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/poly1d.hpp b/models/poly1d.hpp
new file mode 100644
index 0000000..a8c994b
--- /dev/null
+++ b/models/poly1d.hpp
@@ -0,0 +1,67 @@
+#ifndef POLY_MODEL_H_
+#define POLY_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <sstream>
+#include <cassert>
+
+namespace opt_utilities
+{
+ template <typename T,int n>
+ class poly1d
+ :public model<T,T,std::vector<T>,std::string>
+ {
+ private:
+ model<T,T,std::vector<T> >* do_clone()const
+ {
+ return new poly1d<T,n>(*this);
+ }
+ public:
+ poly1d()
+ {
+ assert(n>=0);
+ for(int i=0;i<=n;++i)
+ {
+ std::ostringstream ostr;
+ ostr<<"a"<<i;
+ this->push_param_info(param_info<std::vector<T> >(ostr.str().c_str(),1));
+ }
+
+ }
+
+ public:
+ T do_eval(const T& x,const std::vector<T>& param)
+ {
+ // return x*get_element(param,0)+get_element(param,1);
+ T result(0);
+ for(int i=0;i<=n;++i)
+ {
+ T xn(1);
+ for(int j=0;j<i;++j)
+ {
+ xn*=x;
+ }
+ result+=get_element(param,i)*xn;
+ }
+ return result;
+ }
+
+ private:
+ std::string do_to_string()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
diff --git a/models/pow_model.hpp b/models/pow_model.hpp
new file mode 100644
index 0000000..79baa4b
--- /dev/null
+++ b/models/pow_model.hpp
@@ -0,0 +1,120 @@
+#ifndef POW_MODEL_H_
+#define POW_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename Ty,typename Tx,typename Tp,typename Tstr>
+ class pow_model
+ :public model<Ty,Tx,Tp,Tstr>
+ {
+ private:
+ model<Ty,Tx,Tp,Tstr>* do_clone()const
+ {
+ return new pow_model<Ty,Tx,Tp,Tstr>(*this);
+ }
+ private:
+ pow_model()
+ {
+ }
+
+ private:
+ model<Ty,Tx,Tp,Tstr>* pm1;
+ typename value_type_trait<Tp>::value_type idx;
+
+ public:
+ pow_model(const model<Ty,Tx,Tp,Tstr>& m1,
+ const typename value_type_trait<Tp>::value_type& index)
+ :pm1(m1.clone()),idx(index)
+ {
+ int np1=m1.get_num_params();
+
+ for(int i=0;i<np1;++i)
+ {
+ param_info<Tp,Tstr> p(m1.get_param_info(i));
+ //param_info<Tp,Tstr> p1(p.get_name(),p.get_default_value());
+ this->push_param_info(p);
+ }
+ }
+
+ pow_model(const pow_model& rhs)
+ :pm1(NULL),idx(0)
+ {
+ int np1(0);
+ if(rhs.pm1)
+ {
+ pm1=rhs.pm1->clone();
+ np1=rhs.pm1->get_num_params();
+ for(int i=0;i<np1;++i)
+ {
+ param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+ param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+ this->push_param_info(p1);
+ }
+ }
+ idx=rhs.idx;
+ }
+
+ pow_model& operator=(const pow_model& rhs)
+ {
+ if(this==&rhs)
+ {
+ return *this;
+ }
+ if(!pm1)
+ {
+ //delete pm1;
+ pm1->destroy();
+ }
+
+ int np1(0);
+ if(rhs.pm1)
+ {
+ pm1=rhs.pm1->clone();
+ np1=rhs.pm1->get_num_params();
+ for(int i=0;i<np1;++i)
+ {
+ param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+ // param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+ this->push_param_info(p);
+ }
+ }
+ idx=rhs.idx;
+
+ return *this;
+ }
+
+ ~pow_model()
+ {
+ if(!pm1)
+ {
+ //delete pm1;
+ pm1->destroy();
+ }
+ }
+
+ public:
+ Ty do_eval(const Tx& x,const Tp& param)
+ {
+ if(!pm1)
+ {
+ throw opt_exception("incomplete model!");
+ }
+
+ return std::pow(pm1->eval(x,param),idx);
+ }
+ };
+
+ template <typename Ty,typename Tx,typename Tp,typename Tstr>
+ pow_model<Ty,Tx,Tp,Tstr> pow(const model<Ty,Tx,Tp,Tstr>& m1,
+ const typename value_type_trait<Tp>::value_type& idx)
+ {
+ return pow_model<Ty,Tx,Tp,Tstr>(m1,idx);
+ }
+};
+
+
+
+#endif
+//EOF
diff --git a/models/strmodel1d.cc b/models/strmodel1d.cc
new file mode 100644
index 0000000..df3cd40
--- /dev/null
+++ b/models/strmodel1d.cc
@@ -0,0 +1,82 @@
+#include "strmodel1d.hpp"
+
+using namespace mu;
+using namespace std;
+using namespace opt_utilities;
+
+strmodel1d* strmodel1d::do_clone()const
+{
+ return new strmodel1d(*this);
+}
+
+strmodel1d::strmodel1d()
+{
+ set_buildin_fun();
+}
+
+strmodel1d::strmodel1d(const strmodel1d& rhs)
+ :expr(rhs.expr),
+ par_names(rhs.par_names),
+ var_name(rhs.var_name),
+ par_vec(rhs.par_vec)
+{
+ set_buildin_fun();
+ set_expr(expr,par_names,var_name);
+}
+
+strmodel1d& strmodel1d::operator=(const strmodel1d& rhs)
+{
+ set_buildin_fun();
+ expr=rhs.expr;
+ par_names=rhs.par_names;
+ var_name=rhs.var_name;
+ par_vec=rhs.par_vec;
+ set_expr(expr,par_names,var_name);
+ return *this;
+}
+
+void strmodel1d::set_buildin_fun()
+{
+ mp.DefineFun("sin",sin);
+ mp.DefineFun("cos",cos);
+ mp.DefineFun("exp",exp);
+ mp.DefineFun("tan",tan);
+ mp.DefineFun("log",log);
+}
+
+
+void strmodel1d::set_expr(const string& _expr,
+ const std::vector<std::string>& _par_names,
+ const std::string& _var_name)
+{
+ expr=_expr;
+ par_names=_par_names;
+ var_name=_var_name;
+ this->clear_param_info();
+ par_vec.resize(par_names.size());
+ mp.ClearVar();
+ // mp.ClearFun();
+ for(int i=0;i<par_vec.size();++i)
+ {
+ mp.DefineVar(par_names[i].c_str(),&par_vec[i]);
+ this->push_param_info(param_info<std::vector<double> >(par_names[i],0));
+ }
+ mp.DefineVar(var_name.c_str(),&x);
+ mp.SetExpr(expr.c_str());
+}
+
+
+double strmodel1d::do_eval(const double& _x,const vector<double>& p)
+{
+ for(int i=0;i<par_vec.size();++i)
+ {
+ //get_element(par_vec,i)=get_element(p,i);
+ set_element(par_vec,i,get_element(p,i));
+ // cout<<par_vec[i]<<" ";
+ }
+ x=_x;
+ // cout<<x<<" ";
+ double result(mp.Eval());
+ // cout<<result<<endl;
+ return result;
+}
diff --git a/models/strmodel1d.hpp b/models/strmodel1d.hpp
new file mode 100644
index 0000000..8d992a0
--- /dev/null
+++ b/models/strmodel1d.hpp
@@ -0,0 +1,39 @@
+#ifndef STRMODEL1D_HPP
+#define STRMODEL1D_HPP
+#include <core/fitter.hpp>
+#include <cmath>
+#include <sstream>
+#include <cassert>
+#include <muparser/muParser.h>
+#include <vector>
+#include <string>
+
+class strmodel1d
+ :public opt_utilities::model<double,double,std::vector<double>,std::string>
+{
+private:
+ mu::Parser mp;
+ strmodel1d* do_clone()const;
+ std::vector<double> par_vec;
+ std::vector<std::string> par_names;
+ std::string expr;
+ std::string var_name;
+ double x;
+ void set_buildin_fun();
+public:
+ double do_eval(const double& x,const std::vector<double>& p);
+ strmodel1d();
+ strmodel1d(const strmodel1d& rhs);
+ strmodel1d& operator=(const strmodel1d& rhs);
+
+
+ void set_expr(const std::string& _expr,
+ const std::vector<std::string>& _par_names,
+ const std::string& _var_name);
+
+};
+
+
+
+#endif
+//EOF
diff --git a/models/vecn.hpp b/models/vecn.hpp
new file mode 100644
index 0000000..cac1436
--- /dev/null
+++ b/models/vecn.hpp
@@ -0,0 +1,62 @@
+#ifndef VECN_HPP
+#define VECN_HPP
+#include <iostream>
+namespace opt_utilities
+{
+ template <typename T,int n>
+ class vecn
+ {
+ public:
+ T data[n];
+ public:
+ T& operator[](int i)
+ {
+ return data[i];
+ }
+
+ const T& operator[](int i)const
+ {
+ return data[i];
+ }
+
+ vecn()
+ {
+ for(int i=0;i<n;++i)
+ {
+ data[i]=0;
+ }
+ }
+
+ };
+
+ template <typename T,int n>
+ std::istream& operator>>(std::istream& is,vecn<T,n>& p)
+ {
+ for(int i=0;i<n;++i)
+ {
+ is>>p[i];
+ // std::cout<<i<<std::endl;
+ }
+ return is;
+ }
+
+
+ template <typename T,int n>
+ std::ostream& operator<<(std::ostream& os,const vecn<T,n>& p)
+ {
+ os<<'[';
+ for(int i=0;i<n;++i)
+ {
+ os<<p[i]<<",";
+ // std::cout<<i<<std::endl;
+ }
+ os<<']';
+ return os;
+ }
+
+
+
+};
+
+#endif
+//EOF
Copyright © 2017-2018 Aaron LI. All Rights Reserved.