diff options
| -rw-r--r-- | math/vector_operation.hpp | 19 | ||||
| -rw-r--r-- | statistics/cstat.hpp | 5 | ||||
| -rw-r--r-- | vmodels/normed_dgauss1d.hpp | 75 | ||||
| -rw-r--r-- | vmodels/normed_gauss1d.hpp | 58 |
4 files changed, 156 insertions, 1 deletions
diff --git a/math/vector_operation.hpp b/math/vector_operation.hpp index 3f5532c..a856847 100644 --- a/math/vector_operation.hpp +++ b/math/vector_operation.hpp @@ -15,6 +15,25 @@ namespace opt_utilities } return result; } + + template <typename T> + T contract(const T& x1,const T& x2,...) + { + return x1*x2; + } + + template <typename pT> + typename element_type_trait<pT>::element_type + contract(const pT& v1,const pT& v2,const typename element_type_trait<pT>::element_type&) + { + typename element_type_trait<pT>::element_type result(0); + for(int i=0;i<get_size(v1);++i) + { + result+=get_element(v1,i)*get_element(v2,i); + } + return result; + } + } #endif diff --git a/statistics/cstat.hpp b/statistics/cstat.hpp index 4a647e1..7193a93 100644 --- a/statistics/cstat.hpp +++ b/statistics/cstat.hpp @@ -6,7 +6,10 @@ #define CSTAT_HPP
#define OPT_HEADER
#include <core/fitter.hpp>
+#include <math/vector_operation.hpp>
#include <iostream>
+#include <cassert>
+
using std::cout;using std::endl;
namespace opt_utilities
@@ -50,7 +53,7 @@ namespace opt_utilities for(int i=(this->get_data_set()).size()-1;i>=0;--i)
{
Ty model_y=eval_model(this->get_data_set().get_data(i).get_x(),p);
- result-=this->get_data_set().get_data(i).get_y()*std::log(model_y);
+ result-=contract(this->get_data_set().get_data(i).get_y(),std::log(model_y),result);
}
if(verb)
diff --git a/vmodels/normed_dgauss1d.hpp b/vmodels/normed_dgauss1d.hpp new file mode 100644 index 0000000..4060fa6 --- /dev/null +++ b/vmodels/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/vmodels/normed_gauss1d.hpp b/vmodels/normed_gauss1d.hpp new file mode 100644 index 0000000..c95b324 --- /dev/null +++ b/vmodels/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 |