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 --- mass_profile/Makefile | 86 ------ mass_profile/README.md | 37 --- mass_profile/beta.hpp | 45 --- mass_profile/beta_cfg.cpp | 85 ------ mass_profile/beta_cfg.hpp | 23 -- mass_profile/calc_lx_beta.cpp | 497 ------------------------------ mass_profile/calc_lx_dbeta.cpp | 548 --------------------------------- mass_profile/chisq.hpp | 206 ------------- mass_profile/dbeta.hpp | 108 ------- mass_profile/dump_fit_qdp.cpp | 55 ---- mass_profile/dump_fit_qdp.hpp | 16 - mass_profile/fit_beta_sbp.cpp | 534 -------------------------------- mass_profile/fit_dbeta_sbp.cpp | 586 ------------------------------------ mass_profile/fit_nfw_mass.cpp | 159 ---------- mass_profile/fit_wang2012_model.cpp | 195 ------------ mass_profile/nfw.hpp | 56 ---- mass_profile/projector.hpp | 214 ------------- mass_profile/report_error.cpp | 39 --- mass_profile/report_error.hpp | 6 - mass_profile/spline.hpp | 110 ------- mass_profile/vchisq.hpp | 137 --------- mass_profile/wang2012_model.hpp | 67 ----- 22 files changed, 3809 deletions(-) delete mode 100644 mass_profile/Makefile delete mode 100644 mass_profile/README.md delete mode 100644 mass_profile/beta.hpp delete mode 100644 mass_profile/beta_cfg.cpp delete mode 100644 mass_profile/beta_cfg.hpp delete mode 100644 mass_profile/calc_lx_beta.cpp delete mode 100644 mass_profile/calc_lx_dbeta.cpp delete mode 100644 mass_profile/chisq.hpp delete mode 100644 mass_profile/dbeta.hpp delete mode 100644 mass_profile/dump_fit_qdp.cpp delete mode 100644 mass_profile/dump_fit_qdp.hpp delete mode 100644 mass_profile/fit_beta_sbp.cpp delete mode 100644 mass_profile/fit_dbeta_sbp.cpp delete mode 100644 mass_profile/fit_nfw_mass.cpp delete mode 100644 mass_profile/fit_wang2012_model.cpp delete mode 100644 mass_profile/nfw.hpp delete mode 100644 mass_profile/projector.hpp delete mode 100644 mass_profile/report_error.cpp delete mode 100644 mass_profile/report_error.hpp delete mode 100644 mass_profile/spline.hpp delete mode 100644 mass_profile/vchisq.hpp delete mode 100644 mass_profile/wang2012_model.hpp (limited to 'mass_profile') diff --git a/mass_profile/Makefile b/mass_profile/Makefile deleted file mode 100644 index 17b7a0e..0000000 --- a/mass_profile/Makefile +++ /dev/null @@ -1,86 +0,0 @@ -# -# 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) diff --git a/mass_profile/README.md b/mass_profile/README.md deleted file mode 100644 index 0c61821..0000000 --- a/mass_profile/README.md +++ /dev/null @@ -1,37 +0,0 @@ -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/mass_profile/beta.hpp b/mass_profile/beta.hpp deleted file mode 100644 index 6ae70ac..0000000 --- a/mass_profile/beta.hpp +++ /dev/null @@ -1,45 +0,0 @@ -#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/mass_profile/beta_cfg.hpp b/mass_profile/beta_cfg.hpp deleted file mode 100644 index 27148df..0000000 --- a/mass_profile/beta_cfg.hpp +++ /dev/null @@ -1,23 +0,0 @@ -#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/mass_profile/calc_lx_beta.cpp b/mass_profile/calc_lx_beta.cpp deleted file mode 100644 index 11d7b04..0000000 --- a/mass_profile/calc_lx_beta.cpp +++ /dev/null @@ -1,497 +0,0 @@ -/** - * 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/mass_profile/fit_beta_sbp.cpp b/mass_profile/fit_beta_sbp.cpp deleted file mode 100644 index 295fa1e..0000000 --- a/mass_profile/fit_beta_sbp.cpp +++ /dev/null @@ -1,534 +0,0 @@ -/* - 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/mass_profile/projector.hpp b/mass_profile/projector.hpp deleted file mode 100644 index 4ef4b32..0000000 --- a/mass_profile/projector.hpp +++ /dev/null @@ -1,214 +0,0 @@ -#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/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 -#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