diff options
Diffstat (limited to 'mass_profile/vchisq.hpp')
| -rw-r--r-- | mass_profile/vchisq.hpp | 137 |
1 files changed, 0 insertions, 137 deletions
diff --git a/mass_profile/vchisq.hpp b/mass_profile/vchisq.hpp deleted file mode 100644 index 8cd8481..0000000 --- a/mass_profile/vchisq.hpp +++ /dev/null @@ -1,137 +0,0 @@ -/**
- \file vchisq.hpp
- \brief chi-square statistic
- \author Junhua Gu
- */
-
-#ifndef VCHI_SQ_HPP
-#define VCHI_SQ_HPP
-
-#define OPT_HEADER
-
-#include <core/fitter.hpp>
-#include <iostream>
-#include <vector>
-#include <misc/optvec.hpp>
-#include <cmath>
-
-using std::cerr;
-using std::endl;
-
-namespace opt_utilities
-{
- template<typename T>
- class vchisq
- :public statistic<std::vector<T>,std::vector<T>,std::vector<T>,T,std::string>
- {
- private:
- bool verb;
- int n;
- bool limit_bound;
- typedef std::vector<T> Tp;
-
- vchisq<T>* do_clone()const
- {
- return new vchisq<T>(*this);
- }
-
- const char* do_get_type_name()const
- {
- return "chi^2 statistic";
- }
-
- public:
- void verbose(bool v)
- {
- verb=v;
- }
-
- void set_limit()
- {
- limit_bound=true;
- }
-
- void clear_limit()
- {
- limit_bound=false;
- }
- public:
- vchisq()
- :verb(false),limit_bound(false)
- {}
-
- T do_eval(const std::vector<T>& p)
- {
- if(limit_bound)
- {
- if(!this->get_fitter().get_model().meets_constraint(p))
- {
- return 1e99;
- }
- }
- T result(0);
- std::vector<float> vx;
- std::vector<float> vy;
- std::vector<float> vye;
- std::vector<float> my;
- float x1=1e99,x2=-1e99,y1=1e99,y2=-1e99;
- if(verb)
- {
- n++;
- if(n%100==0)
- {
- vx.resize(this->get_data_set().get_data(0).get_y().size());
- vy.resize(this->get_data_set().get_data(0).get_y().size());
- vye.resize(this->get_data_set().get_data(0).get_y().size());
- my.resize(this->get_data_set().get_data(0).get_y().size());
- }
- }
- for(int i=(this->get_data_set()).size()-1;i>=0;--i)
- {
- const std::vector<double> y_model(this->eval_model(this->get_data_set().get_data(i).get_x(),p));
- const std::vector<double>& y=this->get_data_set().get_data(i).get_y();
- const std::vector<double>& ye=this->get_data_set().get_data(i).get_y_lower_err();
- for(size_t j=0;j<y.size();++j)
- {
- double chi=(y_model[j]-y[j])/ye[j];
- result+=chi*chi;
- }
-
- if(verb&&n%100==0)
- {
- for(size_t j=0;j<y.size();++j)
- {
- vx.at(j)=((this->get_data_set().get_data(i).get_x().at(j)+this->get_data_set().get_data(i).get_x().at(j+1))/2.);
- vy.at(j)=(y[j]);
- vye.at(j)=ye[j];
- my.at(j)=(y_model[j]);
- x1=std::min(vx.at(j),x1);
- y1=std::min(vy.at(j),y1);
- x2=std::max(vx.at(j),x2);
- y2=std::max(vy.at(j),y2);
- vye[j]=log10(vy[j]+vye[j])-log10(vy[j]);
- vx[j]=log10(vx[j]);
- vy[j]=log10(vy[j]);
- my[j]=log10(my[j]);
- }
- }
- }
- if(verb)
- {
- if(n%100==0)
- {
- cerr<<result<<"\t";
- for(size_t i=0;i<get_size(p);++i)
- {
- cerr<<get_element(p,i)<<",";
- }
- cerr<<endl;
- }
- }
-
- return result;
- }
- };
-
-}
-#endif
|