diff options
Diffstat (limited to 'mass_profile/fit_wang2012_model.cpp')
| -rw-r--r-- | mass_profile/fit_wang2012_model.cpp | 195 |
1 files changed, 195 insertions, 0 deletions
diff --git a/mass_profile/fit_wang2012_model.cpp b/mass_profile/fit_wang2012_model.cpp new file mode 100644 index 0000000..3fe14b6 --- /dev/null +++ b/mass_profile/fit_wang2012_model.cpp @@ -0,0 +1,195 @@ +/* + 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(int 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(int 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; + } + } +} |