From 4ea7a05ea9a7352602f1f48a860fd81c36e8bc04 Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Mon, 20 Feb 2017 12:26:17 +0800 Subject: Rename mass_profile to src; Add install & uninstall to Makefile --- src/Makefile | 98 ++++++++ src/README.md | 37 +++ src/beta.hpp | 45 ++++ src/beta_cfg.cpp | 85 +++++++ src/beta_cfg.hpp | 23 ++ src/calc_lx_beta.cpp | 497 ++++++++++++++++++++++++++++++++++++++ src/calc_lx_dbeta.cpp | 548 ++++++++++++++++++++++++++++++++++++++++++ src/chisq.hpp | 206 ++++++++++++++++ src/dbeta.hpp | 108 +++++++++ src/dump_fit_qdp.cpp | 55 +++++ src/dump_fit_qdp.hpp | 16 ++ src/fit_beta_sbp.cpp | 534 +++++++++++++++++++++++++++++++++++++++++ src/fit_dbeta_sbp.cpp | 586 +++++++++++++++++++++++++++++++++++++++++++++ src/fit_nfw_mass.cpp | 159 ++++++++++++ src/fit_wang2012_model.cpp | 195 +++++++++++++++ src/nfw.hpp | 56 +++++ src/projector.hpp | 214 +++++++++++++++++ src/report_error.cpp | 39 +++ src/report_error.hpp | 6 + src/spline.hpp | 110 +++++++++ src/vchisq.hpp | 137 +++++++++++ src/wang2012_model.hpp | 67 ++++++ 22 files changed, 3821 insertions(+) create mode 100644 src/Makefile create mode 100644 src/README.md create mode 100644 src/beta.hpp create mode 100644 src/beta_cfg.cpp create mode 100644 src/beta_cfg.hpp create mode 100644 src/calc_lx_beta.cpp create mode 100644 src/calc_lx_dbeta.cpp create mode 100644 src/chisq.hpp create mode 100644 src/dbeta.hpp create mode 100644 src/dump_fit_qdp.cpp create mode 100644 src/dump_fit_qdp.hpp create mode 100644 src/fit_beta_sbp.cpp create mode 100644 src/fit_dbeta_sbp.cpp create mode 100644 src/fit_nfw_mass.cpp create mode 100644 src/fit_wang2012_model.cpp create mode 100644 src/nfw.hpp create mode 100644 src/projector.hpp create mode 100644 src/report_error.cpp create mode 100644 src/report_error.hpp create mode 100644 src/spline.hpp create mode 100644 src/vchisq.hpp create mode 100644 src/wang2012_model.hpp (limited to 'src') diff --git a/src/Makefile b/src/Makefile new file mode 100644 index 0000000..76ddeaa --- /dev/null +++ b/src/Makefile @@ -0,0 +1,98 @@ +# +# Makefile for `chandra-acis-analysis/mass_profile` tools +# +# Junhua GU +# Weitian LI +# 2016-06-07 +# + + +CXX ?= g++ +CXXFLAGS += -Wall +CXXFLAGS += -Werror + +ifdef OPENMP + CXXFLAGS += -fopenmp +endif + +ifdef DEBUG + CXXFLAGS += -g +else + CXXFLAGS += -O2 +endif + +OPT_UTIL_INC ?= -I../opt_utilities + +TARGETS= fit_dbeta_sbp fit_beta_sbp fit_wang2012_model \ + fit_nfw_mass calc_lx_dbeta calc_lx_beta +HEADERS= projector.hpp spline.hpp vchisq.hpp + +all: $(TARGETS) + +# NOTE: +# Object/source files should placed *before* libraries (order matters) + +fit_dbeta_sbp: fit_dbeta_sbp.o beta_cfg.o report_error.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + +fit_beta_sbp: fit_beta_sbp.o beta_cfg.o dump_fit_qdp.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + +fit_wang2012_model: fit_wang2012_model.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + +fit_nfw_mass: fit_nfw_mass.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + +calc_lx_dbeta: calc_lx_dbeta.o beta_cfg.o report_error.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + +calc_lx_beta: calc_lx_beta.o beta_cfg.o dump_fit_qdp.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + + +fit_dbeta_sbp.o: fit_dbeta_sbp.cpp $(HEADERS) + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +fit_beta_sbp.o: fit_beta_sbp.cpp beta.hpp $(HEADERS) + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +fit_wang2012_model.o: fit_wang2012_model.cpp wang2012_model.hpp chisq.hpp + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +fit_nfw_mass.o: fit_nfw_mass.cpp nfw.hpp chisq.hpp + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +calc_lx_dbeta.o: calc_lx_dbeta.cpp $(HEADERS) + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +calc_lx_beta.o: calc_lx_beta.cpp beta.hpp $(HEADERS) + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +beta_cfg.o: beta_cfg.cpp beta_cfg.hpp + $(CXX) $(CXXFLAGS) -c $< + +report_error.o: report_error.cpp report_error.hpp + $(CXX) $(CXXFLAGS) -c $< + +dump_fit_qdp.o: dump_fit_qdp.cpp dump_fit_qdp.hpp + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +%.o: %.cpp + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + + +clean: + rm -f *.o $(TARGETS) + + +install: $(TARGETS) + @for f in $(TARGETS); do \ + (cd ../bin && ln -svf ../src/$$f . ); \ + done + + +uninstall: + @for f in $(TARGETS); do \ + rm -fv ../bin/$$f; \ + done diff --git a/src/README.md b/src/README.md new file mode 100644 index 0000000..0c61821 --- /dev/null +++ b/src/README.md @@ -0,0 +1,37 @@ +Mass/Luminosity/Flux Calculation Tools +============================== + +Junhua GU, Weitian LI, Zhenghao ZHU + + +Introduction +------------ +This directory contains the tools/scripts used to fit the temperature +profile and surface brightness file, to calculate the gas density profile +and gravitational mass profile, to calculate luminosity and flux and +other related quantities. + + +NOTE +---- +* Mass calculation references: + + Walker et al. 2012, MNRAS, 422, 3503-3515, + [ADS:2012MNRAS.422.3503W](http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W) + + Ettori et al. 2013, Space Science Reviews, 177, 119-154, + [ADS:2013SSRv..177..119E](http://adsabs.harvard.edu/abs/2013SSRv..177..119E) +* Errors are calculated at 68% confident level. +* NFW profile is employed to extrapolate the mass profile. +* We use a self-proposed model (see ``wang2012_model.hpp``) to describe/fit + the gas temperature profile. +* The single-beta and double-beta models (with a constant background) are used + to describe the gas density profile. + + +TODO +---- +* Merge the scripts/tools for single-beta and double-beta SBP models + (to reduce code duplications) +* Uncertainties/errors calculation **maybe** inappropriate/problematic + (e.g., ``analyze_mass_profile.py``); + why not just use the quantile-based (e.g., Q84.15-Q15.85 ~= 68.3%) + uncertainties or standard deviation??? diff --git a/src/beta.hpp b/src/beta.hpp new file mode 100644 index 0000000..6ae70ac --- /dev/null +++ b/src/beta.hpp @@ -0,0 +1,45 @@ +#ifndef BETA +#define BETA +#include "projector.hpp" + +namespace opt_utilities +{ + template + class beta + :public model,std::vector,std::vector > + { + public: + beta() + { + this->push_param_info(param_info,std::string>("n0",1,0,1E99)); + this->push_param_info(param_info,std::string>("beta",.66,0,1E99)); + this->push_param_info(param_info,std::string>("rc",100,0,1E99)); + } + + public: + beta* do_clone()const + { + return new beta(*this); + } + + std::vector do_eval(const std::vector & x, + const std::vector& p) + { + T n0=std::abs(p[0]); + T beta=p[1]; + T rc=p[2]; + + std::vector result(x.size()-1); + for(size_t i=1;i +#include +using namespace std; + +cfg_map parse_cfg_file(std::istream& is) +{ + cfg_map result; + result.rmin_pixel=-1; + result.rmin_kpc=-1; + for(;;) + { + std::string line; + getline(is,line); + line+="\n"; + if(!is.good()) + { + break; + } + string key; + istringstream iss(line); + iss>>key; + if(key=="sbp_data") + { + string value; + iss>>value; + result.sbp_data=value; + } + else if(key=="cfunc_profile") + { + string value; + iss>>value; + result.cfunc_profile=value; + } + else if(key=="tprofile") + { + string value; + iss>>value; + result.tprofile=value; + } + else if(key=="z") + { + double z; + iss>>z; + result.z=z; + } + else if(key=="cm_per_pixel") + { + double cm_per_pixel; + iss>>cm_per_pixel; + result.cm_per_pixel=cm_per_pixel; + } + else if(key=="rmin_pixel") + { + double v; + iss>>v; + result.rmin_pixel=v; + } + else if(key=="rmin_kpc") + { + double v; + iss>>v; + result.rmin_kpc=v; + } + else + { + std::vector value; + for(;;) + { + double v; + iss>>v; + if(!iss.good()) + { + break; + } + value.push_back(v); + } + if(!value.empty()) + { + result.param_map[key]=value; + } + } + } + return result; +} diff --git a/src/beta_cfg.hpp b/src/beta_cfg.hpp new file mode 100644 index 0000000..27148df --- /dev/null +++ b/src/beta_cfg.hpp @@ -0,0 +1,23 @@ +#ifndef BETA_CFG +#define BETA_CFG + +#include +#include +#include +#include + +struct cfg_map +{ + std::string sbp_data; + std::string cfunc_profile; + std::string tprofile; + double z; + double cm_per_pixel; + double rmin_kpc; + double rmin_pixel; + std::map > param_map; +}; + +cfg_map parse_cfg_file(std::istream& is); + +#endif diff --git a/src/calc_lx_beta.cpp b/src/calc_lx_beta.cpp new file mode 100644 index 0000000..11d7b04 --- /dev/null +++ b/src/calc_lx_beta.cpp @@ -0,0 +1,497 @@ +/** + * Calculate the total luminosity and flux within the specified radius. + * + * Base on 'fit_beta_sbp.cpp' and supersede 'calc_lx.cpp' + * + * Author: Junhua Gu + */ + +#include +#include +#include +#include +#include +#include "beta_cfg.hpp" +#include "dump_fit_qdp.hpp" +#include "report_error.hpp" +#include "vchisq.hpp" +#include "chisq.hpp" +#include "beta.hpp" +#include "models/beta1d.hpp" +#include +#include +#include +#include +#include "spline.hpp" + +using namespace std; +using namespace opt_utilities; +//double s=5.63136645E20; +const double kpc=3.086E21;//kpc in cm +const double Mpc=kpc*1000; + +double beta_func(double r, double n0, double rc, double beta) +{ + return abs(n0) * pow(1+r*r/rc/rc, -3./2.*abs(beta)); +} + +//A class enclosing the spline interpolation method +class spline_func_obj + :public func_obj +{ + //has an spline object + spline spl; +public: + //This function is used to calculate the intepolated value + double do_eval(const double& x) + { + return spl.get_value(x); + } + + //we need this function, when this object is performing a clone of itself + spline_func_obj* do_clone()const + { + return new spline_func_obj(*this); + } + +public: + //add points to the spline object, after which the spline will be initialized + void add_point(double x,double y) + { + spl.push_point(x,y); + } + + //before getting the intepolated value, the spline should be initialzied by calling this function + void gen_spline() + { + spl.gen_spline(0,0); + } +}; + +int main(int argc,char* argv[]) +{ + if(argc<4) + { + cerr< [cfunc2_erg ...]"< radii; + std::vector sbps; + std::vector sbpe; + std::vector radii_all; + std::vector sbps_all; + std::vector sbpe_all; + //prepend the zero point to radius list + radii.push_back(0.0); + radii_all.push_back(0.0); + //read sbp and sbp error data + ifstream ifs(cfg.sbp_data.c_str()); + std::string line; + if (ifs.is_open()) + { + while(std::getline(ifs, line)) + { + if (line.empty()) + continue; + + std::istringstream iss(line); + double x, xe, y, ye; + if ((iss >> x >> xe >> y >> ye)) + { + std::cerr << "sbprofile data: " + << x << ", " << xe << ", " << y << ", " << ye + << std::endl; + radii.push_back(x+xe); /* NOTE: use outer radii of regions */ + radii_all.push_back(x+xe); + sbps.push_back(y); + sbps_all.push_back(y); + sbpe.push_back(ye); + sbpe_all.push_back(ye); + } + else + { + std::cerr << "skipped line: " << line << std::endl; + } + } + } + else + { + std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str() + << std::endl; + return 1; + } + + //initialize the cm/pixel value + double cm_per_pixel=cfg.cm_per_pixel; + double rmin; + if(cfg.rmin_pixel>0) + { + rmin=cfg.rmin_pixel; + } + else + { + rmin=cfg.rmin_kpc*kpc/cm_per_pixel; + } + + cerr<<"rmin="< radii_tmp,sbps_tmp,sbpe_tmp; + radii_tmp.resize(radii.size()); + sbps_tmp.resize(sbps.size()); + sbpe_tmp.resize(sbpe.size()); + copy(radii.begin(),radii.end(),radii_tmp.begin()); + copy(sbps.begin(),sbps.end(),sbps_tmp.begin()); + copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin()); + for(list::iterator i=radii_tmp.begin();i!=radii_tmp.end();) + { + if(*i>x>>y; + if(!ifs.good()) + { + break; + } + cerr<radii.back()) + { + break; + } + cf.add_point(x,y); + } + cf.gen_spline(); + + //read temperature profile data + spline_func_obj Tprof; + int tcnt=0; + for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt) + { + assert(ifs1.is_open()); + double x,y; + ifs1>>x>>y; + if(!ifs1.good()) + { + break; + } + cerr<,std::vector > ds; + ds.add_data(data,std::vector >(radii,sbps,sbpe,sbpe,radii,radii)); + + //initial fitter + fitter,vector,vector,double> f; + f.load_data(ds); + //initial the object, which is used to calculate projection effect + projector a; + beta betao; + //attach the cooling function + a.attach_cfunc(cf); + a.set_cm_per_pixel(cm_per_pixel); + a.attach_model(betao); + f.set_model(a); + //chi^2 statistic + vchisq c; + c.verbose(true); + c.set_limit(); + f.set_statistic(c); + //optimization method + f.set_opt_method(powell_method >()); + //initialize the initial values + /* + double n0=0; + double beta=0; + double rc=0; + */ + double bkg_level=0; + + for(std::map >::iterator i=cfg.param_map.begin(); + i!=cfg.param_map.end();++i) + { + std::string pname=i->first; + f.set_param_value(pname,i->second.at(0)); + if(i->second.size()==3) + { + double a1=i->second[1]; + double a2=i->second[2]; + double u=std::max(a1,a2); + double l=std::min(a1,a2); + f.set_param_upper_limit(pname,u); + f.set_param_lower_limit(pname,l); + } + else + { + if(pname=="beta") + { + f.set_param_lower_limit(pname,.3); + f.set_param_upper_limit(pname,1.4); + } + } + } + + f.fit(); + f.fit(); + std::vector p=f.get_all_params(); + /* + n0=f.get_param_value("n0"); + rc=f.get_param_value("rc"); + beta=f.get_param_value("beta"); + */ + //output the datasets and fitting results + ofstream param_output("lx_beta_param.txt"); + for(size_t i=0;i=1) + { + ofs_sbp<<"line on 2"<=1) + { + ofs_sbp<<"no no no"<& pj=dynamic_cast&>(f.get_model()); + pj.attach_cfunc(cf_erg); + + + + mv=f.eval_model_raw(radii,p); + double flux_erg=0; + for(size_t i=0;i +#include +#include +#include +#include "vchisq.hpp" +#include "dbeta.hpp" +#include "beta_cfg.hpp" +#include +#include +#include +#include +#include "spline.hpp" + +using namespace std; +using namespace opt_utilities; +//double s=5.63136645E20; +const double kpc=3.086E21;//kpc in cm +const double Mpc=kpc*1000; + +double dbeta_func(double r, double n01, double rc1, double beta1, + double n02, double rc2, double beta2) +{ + double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1)); + double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2)); + return v1 + v2; +} + + +//A class enclosing the spline interpolation method +class spline_func_obj + :public func_obj +{ + //has an spline object + spline spl; +public: + //This function is used to calculate the intepolated value + double do_eval(const double& x) + { + return spl.get_value(x); + } + + //we need this function, when this object is performing a clone of itself + spline_func_obj* do_clone()const + { + return new spline_func_obj(*this); + } + +public: + //add points to the spline object, after which the spline will be initialized + void add_point(double x,double y) + { + spl.push_point(x,y); + } + + //before getting the intepolated value, the spline should be initialzied by calling this function + void gen_spline() + { + spl.gen_spline(0,0); + } +}; + +int main(int argc,char* argv[]) +{ + if(argc<4) + { + cerr< [cfunc2_erg ...]"< radii; + std::vector sbps; + std::vector sbpe; + std::vector radii_all; + std::vector sbps_all; + std::vector sbpe_all; + //prepend the zero point to radius list + radii.push_back(0.0); + radii_all.push_back(0.0); + //read sbp and sbp error data + ifstream ifs(cfg.sbp_data.c_str()); + std::string line; + if (ifs.is_open()) + { + while(std::getline(ifs, line)) + { + if (line.empty()) + continue; + + std::istringstream iss(line); + double x, xe, y, ye; + if ((iss >> x >> xe >> y >> ye)) + { + std::cerr << "sbprofile data: " + << x << ", " << xe << ", " << y << ", " << ye + << std::endl; + radii.push_back(x+xe); /* NOTE: use outer radii of regions */ + radii_all.push_back(x+xe); + sbps.push_back(y); + sbps_all.push_back(y); + sbpe.push_back(ye); + sbpe_all.push_back(ye); + } + else + { + std::cerr << "skipped line: " << line << std::endl; + } + } + } + else + { + std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str() + << std::endl; + return 1; + } + + //initialize the cm/pixel value + double cm_per_pixel=cfg.cm_per_pixel; + double rmin; + if(cfg.rmin_pixel>0) + { + rmin=cfg.rmin_pixel; + } + else + { + rmin=cfg.rmin_kpc*kpc/cm_per_pixel; + } + + cerr<<"rmin="< radii_tmp,sbps_tmp,sbpe_tmp; + radii_tmp.resize(radii.size()); + sbps_tmp.resize(sbps.size()); + sbpe_tmp.resize(sbpe.size()); + copy(radii.begin(),radii.end(),radii_tmp.begin()); + copy(sbps.begin(),sbps.end(),sbps_tmp.begin()); + copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin()); + for(list::iterator i=radii_tmp.begin();i!=radii_tmp.end();) + { + if(*i>x>>y; + if(!ifs.good()) + { + break; + } + cerr<radii.back()) + { + break; + } + cf.add_point(x,y); + } + cf.gen_spline(); + + //read temperature profile data + spline_func_obj Tprof; + int tcnt=0; + for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt) + { + assert(ifs1.is_open()); + double x,y; + ifs1>>x>>y; + if(!ifs1.good()) + { + break; + } + cerr<,std::vector > ds; + ds.add_data(data,std::vector >(radii,sbps,sbpe,sbpe,radii,radii)); + + //initial fitter + fitter,vector,vector,double> f; + f.load_data(ds); + //initial the object, which is used to calculate projection effect + projector a; + bool tie_beta=false; + if(cfg.param_map.find("beta")!=cfg.param_map.end() + &&cfg.param_map.find("beta1")==cfg.param_map.end() + &&cfg.param_map.find("beta2")==cfg.param_map.end()) + { + dbeta2 dbetao; + a.attach_model(dbetao); + tie_beta=true; + } + else if((cfg.param_map.find("beta1")!=cfg.param_map.end() + ||cfg.param_map.find("beta2")!=cfg.param_map.end()) + &&cfg.param_map.find("beta")==cfg.param_map.end()) + { + dbeta dbetao; + a.attach_model(dbetao); + tie_beta=false; + } + else + { + cerr<<"Error, cannot decide whether to tie beta together or let them vary freely!"< c; + c.verbose(true); + c.set_limit(); + f.set_statistic(c); + //optimization method + f.set_opt_method(powell_method >()); + //initialize the initial values + /* + double =0; + double rc1=0; + double n02=0; + double rc2=0; + double beta=0; + */ + double bkg_level=0; + + if(tie_beta) + { + f.set_param_value("beta",.7); + f.set_param_lower_limit("beta",.3); + f.set_param_upper_limit("beta",1.4); + } + else + { + f.set_param_value("beta1",.7); + f.set_param_lower_limit("beta1",.3); + f.set_param_upper_limit("beta1",1.4); + f.set_param_value("beta2",.7); + f.set_param_lower_limit("beta2",.3); + f.set_param_upper_limit("beta2",1.4); + } + for(std::map >::iterator i=cfg.param_map.begin(); + i!=cfg.param_map.end();++i) + { + std::string pname=i->first; + f.set_param_value(pname,i->second.at(0)); + if(i->second.size()==3) + { + double a1=i->second[1]; + double a2=i->second[2]; + double u=std::max(a1,a2); + double l=std::min(a1,a2); + f.set_param_upper_limit(pname,u); + f.set_param_lower_limit(pname,l); + } + else + { + if(pname=="beta"||pname=="beta1"||pname=="beta2") + { + f.set_param_lower_limit(pname,.3); + f.set_param_upper_limit(pname,1.4); + } + } + } + + + + //perform the fitting, first freeze beta1, beta2, rc1, and rc2 + if(tie_beta) + { + f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("beta")+ + freeze_param,std::vector,std::vector,std::string>("rc1")+ + freeze_param,std::vector,std::vector,std::string>("rc2") + ); + } + else + { + f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("beta1")+ + freeze_param,std::vector,std::vector,std::string>("beta2")+ + freeze_param,std::vector,std::vector,std::string>("rc1")+ + freeze_param,std::vector,std::vector,std::string>("rc2") + ); + } + + f.fit(); + + f.clear_param_modifier(); + + //then perform the fitting, freeze beta1 and beta2 + //f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("beta")); + //f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("bkg")); + f.fit(); + //f.clear_param_modifier(); + + //finally thaw all parameters + f.fit(); + + /* + double beta1=0; + double beta2=0; + + n01=f.get_param_value("n01"); + rc1=f.get_param_value("rc1"); + n02=f.get_param_value("n02"); + rc2=f.get_param_value("rc2"); + if(tie_beta) + { + beta=f.get_param_value("beta"); + beta1=beta; + beta2=beta; + } + else + { + beta1=f.get_param_value("beta1"); + beta2=f.get_param_value("beta2"); + } + */ + + //output the params + ofstream param_output("lx_dbeta_param.txt"); + //output the datasets and fitting results + for(size_t i=0;i& pj=dynamic_cast&>(f.get_model()); + pj.attach_cfunc(cf_erg); + + mv=f.eval_model_raw(radii,p); + double flux_erg=0; + for(size_t i=0;i +#include +#include +#include +#include + +using std::cerr; +using std::endl; + +namespace opt_utilities +{ + static const int display_interval=10; + /** + \brief chi-square statistic + \tparam Ty the return type of model + \tparam Tx the type of the self-var + \tparam Tp the type of model parameter + \tparam Ts the type of the statistic + \tparam Tstr the type of the string used + */ + template + class chisq + :public statistic + { + }; + template<> + class chisq,double,std::string> + :public statistic ,double,std::string> + { + public: + typedef double Ty; + typedef double Tx; + typedef std::vector Tp; + typedef double Ts; + typedef std::string Tstr; + private: + bool verb; + bool limit_bound; + int n; + + statistic* do_clone()const + { + // return const_cast*>(this); + return new chisq(*this); + } + + const char* do_get_type_name()const + { + return "chi^2 statistics (specialized for double)"; + } + public: + void verbose(bool v) + { + verb=v; + } + + void set_limit() + { + limit_bound=true; + } + + void clear_limit() + { + limit_bound=false; + } + public: + chisq() + :verb(true),limit_bound(false) + {} + + Ty do_eval(const Tp& p) + { + if(limit_bound) + { + Tp p1=this->get_fitter().get_model().reform_param(p); + for(size_t i=0;ithis->get_fitter().get_param_info(i).get_upper_limit()|| + p1[i]get_fitter().get_param_info(i).get_lower_limit()) + { + return 1e99; + } + } + } + Ty result(0); + std::vector vx; + std::vector vy; + std::vector vye1; + std::vector vye2; + std::vector my; + float xmin=1e99,xmax=-1e99,ymin=1e99,ymax=-1e99; + if(verb) + { + n++; + if(n%display_interval==0) + { + vx.resize(this->get_data_set().size()); + vy.resize(this->get_data_set().size()); + vye1.resize(this->get_data_set().size()); + vye2.resize(this->get_data_set().size()); + my.resize(this->get_data_set().size()); + } + } + + for(int i=(this->get_data_set()).size()-1;i>=0;--i) + { + +#ifdef HAVE_X_ERROR + Tx x1=this->get_data_set().get_data(i).get_x()-this->get_data_set().get_data(i).get_x_lower_err(); + Tx x2=this->get_data_set().get_data(i).get_x()+this->get_data_set().get_data(i).get_x_upper_err(); + Tx x=this->get_data_set().get_data(i).get_x(); + Ty errx1=(this->eval_model(x1,p)-this->eval_model(x,p)); + Ty errx2=(this->eval_model(x2,p)-this->eval_model(x,p)); + //Ty errx=0; +#else + Ty errx1=0; + Ty errx2=0; +#endif + + Ty y_model=this->eval_model(this->get_data_set().get_data(i).get_x(),p); + Ty y_obs=this->get_data_set().get_data(i).get_y(); + Ty y_err; + + Ty errx=0; + if(errx10?errx1:-errx1; + } + else + { + errx=errx2>0?errx2:-errx2; + } + } + else + { + if(y_obs0?errx2:-errx2; + } + else + { + errx=errx1>0?errx1:-errx1; + } + } + + if(y_model>y_obs) + { + y_err=this->get_data_set().get_data(i).get_y_upper_err(); + } + else + { + y_err=this->get_data_set().get_data(i).get_y_lower_err(); + } + + Ty chi=(y_obs-y_model)/std::sqrt(y_err*y_err+errx*errx); + + result+=chi*chi; + + if(verb&&n%display_interval==0) + { + vx.at(i)=this->get_data_set().get_data(i).get_x(); + vy.at(i)=this->get_data_set().get_data(i).get_y(); + vye1.at(i)=std::abs(this->get_data_set().get_data(i).get_y_lower_err()); + vye2.at(i)=std::abs(this->get_data_set().get_data(i).get_y_upper_err()); + my.at(i)=y_model; + + xmin=std::min(vx.at(i),xmin); + ymin=std::min(vy.at(i),ymin-vye1[i]); + xmax=std::max(vx.at(i),xmax); + ymax=std::max(vy.at(i),ymax+vye2[i]); + } + + } + if(verb) + { + if(n%display_interval==0) + { + cerr< + class dbeta + :public model,std::vector,std::vector > + { + public: + dbeta() + { + this->push_param_info(param_info,std::string>("n01",1)); + this->push_param_info(param_info,std::string>("beta1",.66)); + this->push_param_info(param_info,std::string>("rc1",100)); + + this->push_param_info(param_info,std::string>("n02",1)); + this->push_param_info(param_info,std::string>("beta2",.67)); + this->push_param_info(param_info,std::string>("rc2",110)); + + } + + public: + dbeta* do_clone()const + { + return new dbeta(*this); + } + + std::vector do_eval(const std::vector & x, + const std::vector& p) + { + T n01=std::abs(p[0]); + T beta1=p[1]; + T rc1=p[2]; + + T n02=std::abs(p[3]); + T beta2=p[4]; + T rc2=p[5]; + + + + std::vector result(x.size()-1); + for(size_t i=1;i + class dbeta2 + :public model,std::vector,std::vector > + { + public: + dbeta2() + { + this->push_param_info(param_info,std::string>("n01",1)); + this->push_param_info(param_info,std::string>("rc1",100)); + this->push_param_info(param_info,std::string>("n02",1)); + this->push_param_info(param_info,std::string>("rc2",110)); + this->push_param_info(param_info,std::string>("beta",.67)); + + } + + public: + dbeta2* do_clone()const + { + return new dbeta2(*this); + } + + std::vector do_eval(const std::vector & x, + const std::vector& p) + { + T n01=std::abs(p[0]); + T rc1=p[1]; + + T n02=std::abs(p[2]); + T rc2=p[3]; + T beta=p[4]; + T beta1=beta; + T beta2=beta; + + std::vector result(x.size()-1); + for(size_t i=1;i,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& y,const std::vector& ye) + { + os<<"read serr 1 2"< p=f.get_all_params(); + for(size_t i=1;i,std::vector,std::vector,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& sbps,const std::vector& sbpe) + { + os<<"read serr 1 2"< p=f.get_all_params(); + std::vector mv=f.eval_model_raw(r,p); + for(size_t i=1;i +#include +#include +#include + +namespace opt_utilities +{ + void dump_sbp_beta(std::ostream& os,fitter,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& y,const std::vector& ye); + void dump_rho_beta(std::ostream& os,fitter,std::vector,std::vector,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& sbps,const std::vector& sbpe); + void dump_rho_dbeta(std::ostream& os,fitter,std::vector,std::vector,double,std::string>& f,double cm_per_pixel); +}; + +#endif diff --git a/src/fit_beta_sbp.cpp b/src/fit_beta_sbp.cpp new file mode 100644 index 0000000..295fa1e --- /dev/null +++ b/src/fit_beta_sbp.cpp @@ -0,0 +1,534 @@ +/* + Perform a double-beta density model fitting to the surface brightness data + Author: Junhua Gu + Last modified: 2016.06.07 + This code is distributed with no warrant +*/ + +#include +#include +#include +#include +#include +#include "beta_cfg.hpp" +#include "dump_fit_qdp.hpp" +#include "report_error.hpp" +#include "vchisq.hpp" +#include "chisq.hpp" +#include "beta.hpp" +#include "models/beta1d.hpp" +#include +#include +#include +#include +#include "spline.hpp" + +using namespace std; +using namespace opt_utilities; +//double s=5.63136645E20; +const double kpc=3.086E21;//kpc in cm +const double Mpc=kpc*1000; + +double beta_func(double r, double n0, double rc, double beta) +{ + return abs(n0) * pow(1+r*r/rc/rc, -3./2.*abs(beta)); +} + +//calculate critical density from z, under following cosmological constants +static double calc_critical_density(double z, + const double H0=2.3E-18, + const double Omega_m=.27) +{ + const double G=6.673E-8;//cm^3 g^-1 s^2 + const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m); + const double H=H0*E; + return 3*H*H/8/pi/G; +} + + +//A class enclosing the spline interpolation method +class spline_func_obj + :public func_obj +{ + //has an spline object + spline spl; +public: + //This function is used to calculate the intepolated value + double do_eval(const double& x) + { + return spl.get_value(x); + } + + //we need this function, when this object is performing a clone of itself + spline_func_obj* do_clone()const + { + return new spline_func_obj(*this); + } + +public: + //add points to the spline object, after which the spline will be initialized + void add_point(double x,double y) + { + spl.push_point(x,y); + } + + //before getting the intepolated value, the spline should be initialzied by calling this function + void gen_spline() + { + spl.gen_spline(0,0); + } +}; + +int main(int argc,char* argv[]) +{ + if(argc!=2) + { + cerr<"< radii; + std::vector sbps; + std::vector sbpe; + std::vector radii_all; + std::vector sbps_all; + std::vector sbpe_all; + //prepend the zero point to radius list + radii.push_back(0.0); + radii_all.push_back(0.0); + //read sbp and sbp error data + cerr << "Read surface brightness profile data ..." << endl; + ifstream ifs(cfg.sbp_data.c_str()); + std::string line; + if (ifs.is_open()) + { + while(std::getline(ifs, line)) + { + if (line.empty()) + continue; + + std::istringstream iss(line); + double x, xe, y, ye; + if ((iss >> x >> xe >> y >> ye)) + { + std::cerr << "sbprofile data: " + << x << ", " << xe << ", " << y << ", " << ye + << std::endl; + radii.push_back(x+xe); /* NOTE: use outer radii of regions */ + radii_all.push_back(x+xe); + sbps.push_back(y); + sbps_all.push_back(y); + sbpe.push_back(ye); + sbpe_all.push_back(ye); + } + else + { + std::cerr << "skipped line: " << line << std::endl; + } + } + } + else + { + std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str() + << std::endl; + return 1; + } + + //initialize the cm/pixel value + double cm_per_pixel=cfg.cm_per_pixel; + double rmin; + if(cfg.rmin_pixel>0) + { + rmin=cfg.rmin_pixel; + } + else + { + rmin=cfg.rmin_kpc*kpc/cm_per_pixel; + } + + cerr<<"rmin="< radii_tmp,sbps_tmp,sbpe_tmp; + radii_tmp.resize(radii.size()); + sbps_tmp.resize(sbps.size()); + sbpe_tmp.resize(sbpe.size()); + copy(radii.begin(),radii.end(),radii_tmp.begin()); + copy(sbps.begin(),sbps.end(),sbps_tmp.begin()); + copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin()); + for(list::iterator i=radii_tmp.begin();i!=radii_tmp.end();) + { + if(*i>x>>y; + if(!ifs.good()) + { + break; + } + cerr<radii.back()) + { + cerr << "radius_max: " << radii.back() << endl; + break; + } + cf.add_point(x,y); + } + cf.gen_spline(); + + //read temperature profile data + cerr << "Read temperature profile data ..." << endl; + spline_func_obj Tprof; + int tcnt=0; + for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt) + { + assert(ifs1.is_open()); + double x,y; + ifs1>>x>>y; + if(!ifs1.good()) + { + break; + } + cerr<,std::vector > ds; + ds.add_data(data,std::vector >(radii,sbps,sbpe,sbpe,radii,radii)); + + //initial fitter + fitter,vector,vector,double> f; + f.load_data(ds); + //initial the object, which is used to calculate projection effect + projector a; + beta betao; + //attach the cooling function + a.attach_cfunc(cf); + a.set_cm_per_pixel(cm_per_pixel); + a.attach_model(betao); + f.set_model(a); + //chi^2 statistic + vchisq c; + c.verbose(true); + c.set_limit(); + f.set_statistic(c); + //optimization method + f.set_opt_method(powell_method >()); + //initialize the initial values + double n0=0; + //double beta=atof(arg_map["beta"].c_str()); + double beta=0; + double rc=0; + double bkg_level=0; + + for(std::map >::iterator i=cfg.param_map.begin(); + i!=cfg.param_map.end();++i) + { + std::string pname=i->first; + f.set_param_value(pname,i->second.at(0)); + if(i->second.size()==3) + { + double a1=i->second[1]; + double a2=i->second[2]; + double u=std::max(a1,a2); + double l=std::min(a1,a2); + f.set_param_upper_limit(pname,u); + f.set_param_lower_limit(pname,l); + } + else + { + if(pname=="beta") + { + f.set_param_lower_limit(pname,.3); + f.set_param_upper_limit(pname,1.4); + } + } + } + + f.fit(); + f.fit(); + std::vector p=f.get_all_params(); + n0=f.get_param_value("n0"); + rc=f.get_param_value("rc"); + beta=f.get_param_value("beta"); + //output the datasets and fitting results + ofstream param_output("beta_param.txt"); + for(size_t i=0;i=1) + { + ofs_sbp<<"line on 2"<=1) + { + ofs_sbp<<"no no no"< +#include +#include +#include +#include "vchisq.hpp" +#include "dbeta.hpp" +#include "beta_cfg.hpp" +#include +#include +#include +#include +#include "spline.hpp" + +using namespace std; +using namespace opt_utilities; +//double s=5.63136645E20; +const double kpc=3.086E21;//kpc in cm +const double Mpc=kpc*1000; + +double dbeta_func(double r, double n01, double rc1, double beta1, + double n02, double rc2, double beta2) +{ + double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1)); + double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2)); + return v1 + v2; +} + +//calculate critical density from z, under following cosmological constants +static double calc_critical_density(double z, + const double H0=2.3E-18, + const double Omega_m=.27) +{ + const double G=6.673E-8;//cm^3 g^-1 s^2 + const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m); + const double H=H0*E; + return 3*H*H/8/pi/G; +} + + +//A class enclosing the spline interpolation method +class spline_func_obj + :public func_obj +{ + //has an spline object + spline spl; +public: + //This function is used to calculate the intepolated value + double do_eval(const double& x) + { + return spl.get_value(x); + } + + //we need this function, when this object is performing a clone of itself + spline_func_obj* do_clone()const + { + return new spline_func_obj(*this); + } + +public: + //add points to the spline object, after which the spline will be initialized + void add_point(double x,double y) + { + spl.push_point(x,y); + } + + //before getting the intepolated value, the spline should be initialzied by calling this function + void gen_spline() + { + spl.gen_spline(0,0); + } +}; + +int main(int argc,char* argv[]) +{ + if(argc!=2) + { + cerr<"< radii; + std::vector sbps; + std::vector sbpe; + std::vector radii_all; + std::vector sbps_all; + std::vector sbpe_all; + //prepend the zero point to radius list + radii.push_back(0.0); + radii_all.push_back(0.0); + //read sbp and sbp error data + ifstream ifs(cfg.sbp_data.c_str()); + std::string line; + if (ifs.is_open()) + { + while(std::getline(ifs, line)) + { + if (line.empty()) + continue; + + std::istringstream iss(line); + double x, xe, y, ye; + if ((iss >> x >> xe >> y >> ye)) + { + std::cerr << "sbprofile data: " + << x << ", " << xe << ", " << y << ", " << ye + << std::endl; + radii.push_back(x+xe); /* NOTE: use outer radii of regions */ + radii_all.push_back(x+xe); + sbps.push_back(y); + sbps_all.push_back(y); + sbpe.push_back(ye); + sbpe_all.push_back(ye); + } + else + { + std::cerr << "skipped line: " << line << std::endl; + } + } + } + else + { + std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str() + << std::endl; + return 1; + } + + //initialize the cm/pixel value + double cm_per_pixel=cfg.cm_per_pixel; + double rmin; + if(cfg.rmin_pixel>0) + { + rmin=cfg.rmin_pixel; + } + else + { + rmin=cfg.rmin_kpc*kpc/cm_per_pixel; + } + + cerr<<"rmin="< radii_tmp,sbps_tmp,sbpe_tmp; + radii_tmp.resize(radii.size()); + sbps_tmp.resize(sbps.size()); + sbpe_tmp.resize(sbpe.size()); + copy(radii.begin(),radii.end(),radii_tmp.begin()); + copy(sbps.begin(),sbps.end(),sbps_tmp.begin()); + copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin()); + for(list::iterator i=radii_tmp.begin();i!=radii_tmp.end();) + { + if(*i>x>>y; + if(!ifs.good()) + { + break; + } + cerr<radii.back()) + { + break; + } + cf.add_point(x,y); + } + cf.gen_spline(); + + //read temperature profile data + spline_func_obj Tprof; + int tcnt=0; + for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt) + { + assert(ifs1.is_open()); + double x,y; + ifs1>>x>>y; + if(!ifs1.good()) + { + break; + } + cerr<,std::vector > ds; + ds.add_data(data,std::vector >(radii,sbps,sbpe,sbpe,radii,radii)); + + //initial fitter + fitter,vector,vector,double> f; + f.load_data(ds); + //initial the object, which is used to calculate projection effect + projector a; + bool tie_beta=false; + if(cfg.param_map.find("beta")!=cfg.param_map.end() + &&cfg.param_map.find("beta1")==cfg.param_map.end() + &&cfg.param_map.find("beta2")==cfg.param_map.end()) + { + dbeta2 dbetao; + a.attach_model(dbetao); + tie_beta=true; + } + else if((cfg.param_map.find("beta1")!=cfg.param_map.end() + ||cfg.param_map.find("beta2")!=cfg.param_map.end()) + &&cfg.param_map.find("beta")==cfg.param_map.end()) + { + dbeta dbetao; + a.attach_model(dbetao); + tie_beta=false; + } + else + { + cerr<<"Error, cannot decide whether to tie beta together or let them vary freely!"< c; + c.verbose(true); + c.set_limit(); + f.set_statistic(c); + //optimization method + f.set_opt_method(powell_method >()); + //initialize the initial values + double n01=0; + double rc1=0; + double n02=0; + double rc2=0; + double beta=0; + double bkg_level=0; + if(tie_beta) + { + f.set_param_value("beta",.7); + f.set_param_lower_limit("beta",.3); + f.set_param_upper_limit("beta",1.4); + } + else + { + f.set_param_value("beta1",.7); + f.set_param_lower_limit("beta1",.3); + f.set_param_upper_limit("beta1",1.4); + f.set_param_value("beta2",.7); + f.set_param_lower_limit("beta2",.3); + f.set_param_upper_limit("beta2",1.4); + } + for(std::map >::iterator i=cfg.param_map.begin(); + i!=cfg.param_map.end();++i) + { + std::string pname=i->first; + f.set_param_value(pname,i->second.at(0)); + if(i->second.size()==3) + { + double a1=i->second[1]; + double a2=i->second[2]; + double u=std::max(a1,a2); + double l=std::min(a1,a2); + f.set_param_upper_limit(pname,u); + f.set_param_lower_limit(pname,l); + } + else + { + if(pname=="beta"||pname=="beta1"||pname=="beta2") + { + f.set_param_lower_limit(pname,.3); + f.set_param_upper_limit(pname,1.4); + } + } + } + + + + //perform the fitting, first freeze beta1, beta2, rc1, and rc2 + if(tie_beta) + { + f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("beta")+ + freeze_param,std::vector,std::vector,std::string>("rc1")+ + freeze_param,std::vector,std::vector,std::string>("rc2") + ); + } + else + { + f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("beta1")+ + freeze_param,std::vector,std::vector,std::string>("beta2")+ + freeze_param,std::vector,std::vector,std::string>("rc1")+ + freeze_param,std::vector,std::vector,std::string>("rc2") + ); + } + + f.fit(); + + f.clear_param_modifier(); + + //then perform the fitting, freeze beta1 and beta2 + //f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("beta")); + //f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("bkg")); + f.fit(); + //f.clear_param_modifier(); + + //finally thaw all parameters + f.fit(); + double beta1=0; + double beta2=0; + + n01=f.get_param_value("n01"); + rc1=f.get_param_value("rc1"); + n02=f.get_param_value("n02"); + rc2=f.get_param_value("rc2"); + if(tie_beta) + { + beta=f.get_param_value("beta"); + beta1=beta; + beta2=beta; + } + else + { + beta1=f.get_param_value("beta1"); + beta2=f.get_param_value("beta2"); + } + //output the params + ofstream param_output("dbeta_param.txt"); + //output the datasets and fitting results + for(size_t i=0;i +#include +#include +#include "chisq.hpp" +#include +#include +#include +#include +#include + +using namespace opt_utilities; +using namespace std; +const double cm=1; +const double kpc=3.08568e+21*cm; +const double pi=4*atan(1); +static double calc_critical_density(double z, + const double H0=2.3E-18, + const double Omega_m=.27) +{ + const double G=6.673E-8;//cm^3 g^-1 s^2 + const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m); + const double H=H0*E; + return 3*H*H/8/pi/G; +} + + +int main(int argc,char* argv[]) +{ + if(argc<3) + { + cerr<<"Usage:"< [rmin in kpc]"<=4) + { + rmin_kpc=atof(argv[3]); + } + double z=0; + z=atof(argv[2]); + //define the fitter + fitter,double,std::string> fit; + //define the data set + default_data_set ds; + //open the data file + ifstream ifs(argv[1]); + //cout<<"read serr 2"<>r>>re>>m>>me; + if(!ifs.good()) + { + break; + } + if(r d(r,m,me,me,re,re); + ofs_fit_result< >()); + //use chi^2 statistic + fit.set_statistic(chisq,double,std::string>()); + fit.set_model(nfw()); + //fit.set_param_value("rs",4); + //fit.set_param_value("rho0",100); + fit.fit(); + fit.fit(); + vector p=fit.fit(); + //output parameters + ofstream ofs_param("nfw_param.txt"); + for(size_t i=0;i d=ds.get_data(i); + double x=d.get_x(); + double y=d.get_y(); + double ye=d.get_y_lower_err(); + double ym=fit.eval_model(x,p); + ofs_fit_result< +#include +#include +#include "chisq.hpp" +#include +#include +#include +#include +#include +#include + +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:"< [param file] [cm per pixel]"<=4) + { + cm_per_pixel=atof(argv[3]); + } + + //define the fitter + fitter,double,std::string> fit; + //define the data set + default_data_set ds; + //open the data file + ifstream ifs(argv[1]); + double min_r=1e9; + //cout<<"read serr 2"<0) + { + ofs_fit_result<<"la x radius (kpc)"<>r>>re>>t>>te; + if(!ifs.good()) + { + break; + } + min_r=min(r,min_r); + data d(r,t,te,te,re,re); + //std::cerr<0) + { + ofs_fit_result< >()); + //use chi^2 statistic + chisq,double,std::string> chisq_object; + chisq_object.set_limit(); + fit.set_statistic(chisq_object); + //fit.set_statistic(chisq,double,std::string>()); + fit.set_model(wang2012_model()); + + if(argc>=3&&std::string(argv[2])!="NONE") + { + std::vector 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"<,std::string> fp(freeze_list[0]); + fit.set_param_modifier(fp); + for(size_t i=1;i,std::string>&>(fit.get_param_modifier())+=freeze_param,std::string>(freeze_list[i]); + } + } + } + + for(int i=0;i<100;++i) + { + fit.fit(); + } + vector 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=3&&std::string(argv[2])!="NONE") +#if 0 + { + output_param<0) + { + ofs_fit_result< +#include + +namespace opt_utilities +{ + template + class nfw + :public model,std::string> + { + private: + model >* do_clone()const + { + return new nfw(*this); + } + + const char* do_get_type_name()const + { + return "1d power law"; + } + public: + nfw() + { + this->push_param_info(param_info >("rho0",1,0,1e99)); + this->push_param_info(param_info >("rs",100,0,1e99)); + } + + T do_eval(const T& r,const std::vector& param) + { + T rho0=std::abs(param[0]); + T rs=std::abs(param[1]); + static const T pi=4*std::atan(1); + return 4*pi*rho0*rs*rs*rs*(std::log((r+rs)/rs)-r/(r+rs)); + } + + private: + std::string do_get_information()const + { + return ""; + } + }; +} + + + +#endif +//EOF diff --git a/src/projector.hpp b/src/projector.hpp new file mode 100644 index 0000000..4ef4b32 --- /dev/null +++ b/src/projector.hpp @@ -0,0 +1,214 @@ +#ifndef PROJ_HPP +#define PROJ_HPP +/* + Defining the class that is used to consider the projection effect + Author: Junhua Gu + Last modified: 2011.01.01 +*/ + + +#include +#include +#include + +static const double pi=4*atan(1); +// Ratio of the electron density (n_e) to the proton density (n_p) +// n_gas = n_e + n_p ~= 1.826 n_e => n_e / n_p ~= 1.21 +// Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below +static const double ne_np_ratio = 1.21; + +namespace opt_utilities +{ + //This is used to project a 3-D surface brightness model to 2-D profile + template + class projector + :public model,std::vector,std::vector > + { + private: + //Points to a 3-D model that is to be projected + model,std::vector,std::vector >* pmodel; + func_obj* pcfunc; + T cm_per_pixel; + public: + //default cstr + projector() + :pmodel(NULL_PTR),pcfunc(NULL_PTR),cm_per_pixel(1) + {} + //copy cstr + projector(const projector& rhs) + :cm_per_pixel(rhs.cm_per_pixel) + { + attach_model(*(rhs.pmodel)); + if(rhs.pcfunc) + { + pcfunc=rhs.pcfunc->clone(); + } + else + { + pcfunc=NULL_PTR; + } + } + //assign operator + projector& operator=(const projector& rhs) + { + cm_per_pixel=rhs.cm_per_pixel; + if(pmodel) + { + pmodel->destroy(); + } + if(pcfunc) + { + pcfunc->destroy(); + } + if(rhs.pcfunc) + { + pcfunc=rhs.pcfunc->clone(); + } + if(rhs.pmodel) + { + pmodel=rhs.pmodel->clone(); + } + } + //destr + ~projector() + { + if(pmodel) + { + pmodel->destroy(); + } + if(pcfunc) + { + pcfunc->destroy(); + } + } + //used to clone self + model,std::vector,std::vector >* + do_clone()const + { + return new projector(*this); + } + + public: + void set_cm_per_pixel(const T& x) + { + cm_per_pixel=x; + } + + //attach the model that is to be projected + void attach_model(const model,std::vector,std::vector >& m) + { + this->clear_param_info(); + for(size_t i=0;ipush_param_info(m.get_param_info(i)); + } + this -> push_param_info(param_info, + std::string>("bkg",0,0,1E99)); + pmodel=m.clone(); + pmodel->clear_param_modifier(); + } + + void attach_cfunc(const func_obj& cf) + { + if(pcfunc) + { + pcfunc->destroy(); + } + pcfunc=cf.clone(); + } + + public: + //calc the volume + /* + This is a sphere that is subtracted by a cycline. + /| |\ + / | | \ + | | | | + | | | | + \ | | / + \| |/ + */ + T calc_v_ring(T rsph,T rcyc) + { + if(rcyc& rlist,int nsph,int nrad) + { + if(nsph& p)const + { + std::vector p1(this->reform_param(p)); + for(size_t i=0;i!=p1.size();++i) + { + if(get_element(p1,i)>this->get_param_info(i).get_upper_limit()|| + get_element(p1,i)get_param_info(i).get_lower_limit()) + { + // std::cerr<get_param_info(i).get_name()<<"\t"< p2(p1.size()-1); + for(size_t i=0;imeets_constraint(p2); + } + public: + //Perform the projection + std::vector do_eval(const std::vector& x,const std::vector& p) + { + T bkg=std::abs(p.back()); + //I think following codes are clear enough :). + std::vector unprojected(pmodel->eval(x,p)); + std::vector projected(unprojected.size()); + + for(size_t nrad=0; nrad + +using namespace std; +void report_error(const char* message) +{ + cerr<<"MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM"< +#include +#include +#include +#include + +template +class spline +{ +public: + std::vector x_list; + std::vector y_list; + std::vector y2_list; + +public: + void push_point(T x,T y) + { + if(!x_list.empty()) + { + assert(x>*(x_list.end()-1)); + } + x_list.push_back(x); + y_list.push_back(y); + } + + T get_value(T x) + { + if(x<=x_list[0]) + { + return y_list[0]; + } + if(x>=x_list.back()) + { + return y_list.back(); + } + assert(x_list.size()==y2_list.size()); + assert(x>x_list[0]); + assert(xx) + { + n2--; + } + } + T h=x_list[n2]-x_list[n1]; + double a=(x_list[n2]-x)/h; + double b=(x-x_list[n1])/h; + return a*y_list[n1]+b*y_list[n2]+((a*a*a-a)*y2_list[n1]+ + (b*b*b-b)*y2_list[n2])*(h*h)/6.; + + } + + void gen_spline(T y2_0,T y2_N) + { + int n=x_list.size(); + y2_list.resize(0); + y2_list.resize(x_list.size()); + std::vector u(x_list.size()); + if(std::abs(y2_0)::epsilon()) + { + y2_list[0]=0; + u[0]=0; + } + else + { + y2_list[0]=-.5; + u[0]=(3./(x_list[1]-x_list[0]))*((y_list[1]-y_list[0])/(x_list[1]-x_list[0])-y2_0); + } + for(int i=1;i::epsilon()) + { + qn=un=0; + } + else + { + qn=.5; + un=(3./(x_list[n-1]-x_list[n-2]))*(y2_N-(y_list[n-1]-y_list[n-2])/(x_list[n-1]-x_list[n-2])); + + } + y2_list[n-1]=(un-qn*u[n-2])/(qn*y2_list[n-2]+1.); + for(int i=n-2;i>=0;--i) + { + y2_list[i]=y2_list[i]*y2_list[i+1]+u[i]; + } + } + +}; + +#endif /* SPLINE_HPP */ diff --git a/src/vchisq.hpp b/src/vchisq.hpp new file mode 100644 index 0000000..8cd8481 --- /dev/null +++ b/src/vchisq.hpp @@ -0,0 +1,137 @@ +/** + \file vchisq.hpp + \brief chi-square statistic + \author Junhua Gu + */ + +#ifndef VCHI_SQ_HPP +#define VCHI_SQ_HPP + +#define OPT_HEADER + +#include +#include +#include +#include +#include + +using std::cerr; +using std::endl; + +namespace opt_utilities +{ + template + class vchisq + :public statistic,std::vector,std::vector,T,std::string> + { + private: + bool verb; + int n; + bool limit_bound; + typedef std::vector Tp; + + vchisq* do_clone()const + { + return new vchisq(*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& p) + { + if(limit_bound) + { + if(!this->get_fitter().get_model().meets_constraint(p)) + { + return 1e99; + } + } + T result(0); + std::vector vx; + std::vector vy; + std::vector vye; + std::vector 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 y_model(this->eval_model(this->get_data_set().get_data(i).get_x(),p)); + const std::vector& y=this->get_data_set().get_data(i).get_y(); + const std::vector& ye=this->get_data_set().get_data(i).get_y_lower_err(); + for(size_t j=0;jget_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< +#include + +namespace opt_utilities +{ + template + class wang2012_model + :public model,std::string> + { + private: + model >* do_clone()const + { + return new wang2012_model(*this); + } + + const char* do_get_type_name()const + { + return "1d power law"; + } + public: + wang2012_model() + { + this->push_param_info(param_info >("A",5,0,500)); + this->push_param_info(param_info >("n",1.66,0,10)); + this->push_param_info(param_info >("xi",0.45,0,1)); + this->push_param_info(param_info >("a2",1500,0,1e8)); + this->push_param_info(param_info >("a3",50,0,1e8)); + this->push_param_info(param_info >("beta",0.49,0.1,0.7)); + this->push_param_info(param_info >("T0",0,0,10)); + + } + + T do_eval(const T& x,const std::vector& param) + { + T A=param[0]; + T n=param[1]; + T xi=param[2]; + T a2=param[3]; + T a3=param[4]; + T beta=param[5]; + T T0=param[6]; + return A*(pow(x,n)+xi*a2)/(pow(x,n)+a2)/pow(1+x*x/a3/a3,beta)+T0; + //return A*(pow(x,n)+a1)/(pow(x,n)+1)/pow(1+x*x/a3/a3,beta)+T0; + } + + private: + std::string do_get_information()const + { + return ""; + } + }; +} + + + +#endif +//EOF -- cgit v1.2.2