diff options
Diffstat (limited to 'mass_profile/fit_wang2012_model.cpp')
| -rw-r--r-- | mass_profile/fit_wang2012_model.cpp | 195 |
1 files changed, 0 insertions, 195 deletions
diff --git a/mass_profile/fit_wang2012_model.cpp b/mass_profile/fit_wang2012_model.cpp deleted file mode 100644 index f15a2b1..0000000 --- a/mass_profile/fit_wang2012_model.cpp +++ /dev/null @@ -1,195 +0,0 @@ -/* - Fitting Jy Wang's temperature profile model - Author: Jingying Wang - Last modification 20120819 - -*/ - -#include "wang2012_model.hpp" -#include <core/optimizer.hpp> -#include <core/fitter.hpp> -#include <data_sets/default_data_set.hpp> -#include "chisq.hpp" -#include <methods/powell/powell_method.hpp> -#include <core/freeze_param.hpp> -#include <iostream> -#include <fstream> -#include <vector> -#include <string> - -using namespace opt_utilities; -using namespace std; -const double cm=1; -const double kpc=3.08568e+21*cm; - -int main(int argc,char* argv[]) -{ - if(argc<2) - { - cerr<<"Usage:"<<argv[0]<<" <data file with 4 columns of x, xe, y, ye> [param file] [cm per pixel]"<<endl; - return -1; - } - double cm_per_pixel=-1; - if(argc>=4) - { - cm_per_pixel=atof(argv[3]); - } - - //define the fitter - fitter<double,double,vector<double>,double,std::string> fit; - //define the data set - default_data_set<double,double> ds; - //open the data file - ifstream ifs(argv[1]); - double min_r=1e9; - //cout<<"read serr 2"<<endl; - ofstream ofs_fit_result("fit_result.qdp"); - ofs_fit_result<<"read serr 1 2"<<endl; - ofs_fit_result<<"skip single"<<endl; - if(cm_per_pixel>0) - { - ofs_fit_result<<"la x radius (kpc)"<<endl; - } - else - { - ofs_fit_result<<"la x radius (pixel)"<<endl; - } - ofs_fit_result<<"la y temperature (keV)"<<endl; - for(;;) - { - //read radius, temperature and error - double r,re,t,te; - ifs>>r>>re>>t>>te; - if(!ifs.good()) - { - break; - } - min_r=min(r,min_r); - data<double,double> d(r,t,te,te,re,re); - //std::cerr<<r<<"\t"<<t<<"\t"<<te<<endl; - if(cm_per_pixel>0) - { - ofs_fit_result<<r*cm_per_pixel/kpc<<"\t"<<re*cm_per_pixel/kpc<<"\t"<<t<<"\t"<<te<<endl; - } - else - { - ofs_fit_result<<r<<"\t"<<re<<"\t"<<t<<"\t"<<te<<endl; - } - ds.add_data(d); - } - ofs_fit_result<<"no no no"<<endl; - //load data - fit.load_data(ds); - //define the optimization method - fit.set_opt_method(powell_method<double,vector<double> >()); - //use chi^2 statistic - chisq<double,double,vector<double>,double,std::string> chisq_object; - chisq_object.set_limit(); - fit.set_statistic(chisq_object); - //fit.set_statistic(chisq<double,double,vector<double>,double,std::string>()); - fit.set_model(wang2012_model<double>()); - - if(argc>=3&&std::string(argv[2])!="NONE") - { - std::vector<std::string> freeze_list; - ifstream ifs_param(argv[2]); - assert(ifs_param.is_open()); - for(;;) - { - string pname; - double pvalue; - double lower,upper; - char param_status; - ifs_param>>pname>>pvalue>>lower>>upper>>param_status; - if(!ifs_param.good()) - { - break; - } - if(param_status=='F') - { - freeze_list.push_back(pname); - } - if(pvalue<=lower||pvalue>=upper) - { - cerr<<"Invalid initial value, central value not enclosed by the lower and upper boundaries, adjust automatically"<<endl; - pvalue=std::max(pvalue,lower); - pvalue=std::min(pvalue,upper); - } - fit.set_param_value(pname,pvalue); - fit.set_param_lower_limit(pname,lower); - fit.set_param_upper_limit(pname,upper); - } - if(!freeze_list.empty()) - { - freeze_param<double,double,std::vector<double>,std::string> fp(freeze_list[0]); - fit.set_param_modifier(fp); - for(size_t i=1;i<freeze_list.size();++i) - { - dynamic_cast<freeze_param<double,double,std::vector<double>,std::string>&>(fit.get_param_modifier())+=freeze_param<double,double,std::vector<double>,std::string>(freeze_list[i]); - } - } - } - - for(int i=0;i<100;++i) - { - fit.fit(); - } - vector<double> p=fit.fit(); -#if 0 - ofstream output_param; - if(argc>=3&&std::string(argv[2])!="NONE") - { - output_param.open(argv[2]); - } - else - { - output_param.open("para0.txt"); - } -#endif - //output parameters - for(size_t i=0;i<fit.get_num_params();++i) - { - std::string pname=fit.get_param_info(i).get_name(); - std::string pstatus=fit.get_model().report_param_status(pname); - - if(pstatus==""||pstatus=="thawed") - { - pstatus="T"; - } - else - { - pstatus="F"; - } - cout<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<"\t"<<fit.get_param_info(i).get_lower_limit()<<"\t"<<fit.get_param_info(i).get_upper_limit()<<"\t"<<pstatus<<endl; - //if(argc>=3&&std::string(argv[2])!="NONE") -#if 0 - { - output_param<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<"\t"<<fit.get_param_info(i).get_lower_limit()<<"\t"<<fit.get_param_info(i).get_upper_limit()<<"\t"<<pstatus<<endl; - } -#endif - } - - //cout<<"T0"<<"\t"<<fit.get_param_value("T0")<<endl; - //cout<<"T1"<<"\t"<<fit.get_param_value("T1")<<endl; - //cout<<"xt"<<"\t"<<fit.get_param_value("xt")<<endl; - //cout<<"eta"<<"\t"<<fit.get_param_value("eta")<<endl; - //dump the data for checking - ofstream ofs_model("wang2012_dump.qdp"); - - - min_r=0; - for(double x=min_r;x<3000;x+=10) - { - double model_value=fit.eval_model_raw(x,p); - - ofs_model<<x<<"\t"<<model_value<<endl; - if(cm_per_pixel>0) - { - ofs_fit_result<<x*cm_per_pixel/kpc<<"\t0\t"<<model_value<<"\t0"<<endl; - } - else - { - ofs_fit_result<<x<<"\t0\t"<<model_value<<"\t0"<<endl; - } - } -} |