diff options
| -rw-r--r-- | mass_profile/Makefile | 18 | ||||
| -rw-r--r-- | mass_profile/fit_lt_bpl.cpp | 336 | ||||
| -rw-r--r-- | mass_profile/fit_lt_pl.cpp | 239 | ||||
| -rw-r--r-- | mass_profile/fit_mt_bpl.cpp | 350 | ||||
| -rw-r--r-- | mass_profile/fit_mt_pl.cpp | 265 |
5 files changed, 2 insertions, 1206 deletions
diff --git a/mass_profile/Makefile b/mass_profile/Makefile index 67e9aae..a2067ad 100644 --- a/mass_profile/Makefile +++ b/mass_profile/Makefile @@ -23,10 +23,8 @@ endif OPT_UTIL_INC ?= -I../opt_utilities -TARGETS= fit_dbeta_sbp fit_beta_sbp \ - fit_wang2012_model \ - fit_nfw_mass fit_mt_pl fit_lt_pl fit_mt_bpl \ - fit_lt_bpl calc_lx_dbeta calc_lx_beta +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) @@ -52,18 +50,6 @@ calc_lx_dbeta: calc_lx_dbeta.o beta_cfg.o report_error.o calc_lx_beta: calc_lx_beta.o beta_cfg.o dump_fit_qdp.o $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) -fit_mt_pl: fit_mt_pl.cpp - $(CXX) $(CXXFLAGS) $< -o $@ $(OPT_UTIL_INC) - -fit_lt_pl: fit_lt_pl.cpp - $(CXX) $(CXXFLAGS) $< -o $@ $(OPT_UTIL_INC) - -fit_mt_bpl: fit_mt_bpl.cpp - $(CXX) $(CXXFLAGS) $< -o $@ $(OPT_UTIL_INC) - -fit_lt_bpl: fit_lt_bpl.cpp - $(CXX) $(CXXFLAGS) $< -o $@ $(OPT_UTIL_INC) - fit_dbeta_sbp.o: fit_dbeta_sbp.cpp constrained_dbeta.hpp $(HEADERS) $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) diff --git a/mass_profile/fit_lt_bpl.cpp b/mass_profile/fit_lt_bpl.cpp deleted file mode 100644 index 0c06b22..0000000 --- a/mass_profile/fit_lt_bpl.cpp +++ /dev/null @@ -1,336 +0,0 @@ -/* - Perform a double-beta density model fitting to the surface brightness data - Author: Junhua Gu - Last modified: 2011.01.01 - This code is distributed with no warrant -*/ - -//#define HAVE_X_ERROR -#include <iomanip> -#include <iostream> -#include <sstream> -#include <fstream> -#include <models/bpl1d.hpp> -#include <models/lin1d.hpp> -#include "statistics/chisq.hpp" -#include "statistics/logchisq.hpp" -#include "statistics/leastsq.hpp" -#include <data_sets/default_data_set.hpp> -#include <methods/powell/powell_method.hpp> -#include <methods/gsl_simplex/gsl_simplex.hpp> -#include <core/freeze_param.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; -const double pi=4*atan(1); -double std_norm_rand() -{ - double u=0; - double v=0; - while(u<=0||v<=0) - { - u=rand()/(double)RAND_MAX; - rand(); - v=rand()/(double)RAND_MAX; - } - double x=std::sqrt(-log(u))*cos(2*pi*v); - return x; -} - -double shuffle_data(double xc,double xl,double xu) -{ - double lxc=log(xc); - double lxl=log(xc-xl)-log(xc); - double lxu=log(xc+xu)-log(xc); - - if(std_norm_rand()>0) - { - double result=std::exp(lxc-std::abs(std_norm_rand()*lxl)); - return result; - } - else - { - double result=std::exp(lxc+std::abs(std_norm_rand()*lxu)); - return result; - } -} - -int main(int argc,char* argv[]) -{ - srand(time(0)); - if(argc!=4) - { - cerr<<"Usage:"<<argv[0]<<" <a 5 column file with T -Terr +Terr L Lerr> <initial broken temperature> <T lower limit>"<<endl; - return -1; - } - double T_lower_limit(atof(argv[3])); - double Tb=atof(argv[2]); - assert(Tb>0); - ifstream ifs_data(argv[1]); - default_data_set<double,double> ds; - ofstream ofs_result("l-t_bpl_result.qdp"); - ofs_result<<"read terr 1 2"<<endl; - ofs_result<<"skip single"<<endl; - ofs_result<<"log"<<endl; - //ofs_result<<"li on 2"<<endl; - ofs_result<<"time off"<<endl; - ofs_result<<"la f"<<endl; - ofs_result<<"la x temperature (keV)"<<endl; - ofs_result<<"la y luminosity (10\\u43\\d erg s\\u-1\\d)"<<endl; - double yunit=1E43; - double sxxl=0; - double s1l=0; - double sxl=0; - double syl=0; - double sxyl=0; - - double sxxu=0; - double s1u=0; - double sxu=0; - double syu=0; - double sxyu=0; - - bool is_first_nonono=true; - - for(;;) - { - double T,Tl,Tu; - double L,Ll,Lu; - std::string line; - getline(ifs_data,line); - //ifs_data>>T>>Tl>>Tu>>M>>Ml>>Mu; - if(!ifs_data.good()) - { - break; - } - line+=" "; - istringstream iss(line); - - if(line[0]=='#') - { - if(!is_first_nonono) - { - ofs_result<<"no no no"<<endl; - } - else - { - is_first_nonono=false; - } - continue; - } - iss>>T>>Tl>>Tu>>L>>Ll; - Lu=Ll; - //std::cerr<<L<<"\t"<<Lerr<<endl; - if(!iss.good()) - { - continue; - } - - if(T<T_lower_limit||L<0) - { - continue; - } - if(std::abs(Lu)+std::abs(Ll)<L*.1) - { - double k=L*.1/(std::abs(Lu)+std::abs(Ll)); - Lu*=k; - Ll*=k; - } - if(std::abs(Ll)>std::abs(L)) - { - continue; - } - Tl=std::abs(Tl); - Tu=std::abs(Tu); - Ll=std::abs(Ll); - Lu=std::abs(Lu); - ofs_result<<T<<"\t"<<-std::abs(Tl)<<"\t"<<+std::abs(Tu)<<"\t"<<L/yunit<<"\t"<<-std::abs(Ll/yunit)<<"\t"<<+std::abs(Lu/yunit)<<endl; - double x=(T); - double y=(L); - double xu=Tu; - double xl=Tl; - - double yu=Lu; - double yl=Ll; - if(T>Tb) - { - sxxl+=log(x)*log(x); - sxl+=log(x); - syl+=log(y); - sxyl+=log(y)*log(x); - s1l+=1; - } - else - { - sxxu+=log(x)*log(x); - sxu+=log(x); - syu+=log(y); - sxyu+=log(y)*log(x); - s1u+=1; - } - data<double,double> d(x,y,std::abs(yl),std::abs(yu), - std::abs(xl),std::abs(xu)); - ds.add_data(d); - } - - double Ml=sxxl*s1l-sxl*sxl; - double Mal=sxyl*s1l-syl*sxl; - double Mbl=sxxl*syl-sxl*sxyl; - double k0l=Mal/Ml; - double b0l=Mbl/Ml; - - double Mu=sxxu*s1u-sxu*sxu; - double Mau=sxyu*s1u-syu*sxu; - double Mbu=sxxu*syu-sxu*sxyu; - double k0u=Mau/Mu; - double b0u=Mbu/Mu; - - double gamma0l=k0l; - double gamma0u=k0u; - - double ampl0l=exp(b0l)*pow(Tb,gamma0l); - double ampl0u=exp(b0u)*pow(Tb,gamma0u);; - - - - ofs_result<<"no no no"<<endl; - fitter<double,double,vector<double>,double,std::string> fit; - fit.set_opt_method(powell_method<double,vector<double> >()); - - fit.set_statistic(logchisq<double,double,vector<double>,double,std::string>()); - //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>()); - fit.set_model(bpl1d<double>()); - fit.load_data(ds); - - cerr<<"k0l="<<k0l<<endl; - cerr<<"k0u="<<k0u<<endl; - cerr<<"Ampl0="<<(ampl0l+ampl0u)/2<<endl; - - fit.set_param_value("bpx",Tb); - fit.set_param_value("bpy",(ampl0l+ampl0u)/2); - fit.set_param_value("gamma1",gamma0l); - fit.set_param_value("gamma2",gamma0u); - - fit.fit(); - //fit.set_opt_method(gsl_simplex<double,vector<double> >()); - fit.fit(); - std::vector<double> p=fit.fit(); - Tb=fit.get_param_value("bpx"); - ///std::cout<<"chi="<<fit.get_statistic().eval(p)<<std::endl; - for(double i=.5;i<12;i*=1.01) - { - ofs_result<<i<<"\t0\t0\t"<<fit.eval_model_raw(i,p)/yunit<<"\t0\t0\n"; - } - - - std::vector<double> mean_p(p.size()); - std::vector<double> mean_p2(p.size()); - int cnt=0; - for(int n=0;n<100;++n) - { - ++cnt; - cerr<<"."; - double sxxl=0; - double s1l=0; - double sxl=0; - double syl=0; - double sxyl=0; - - double sxxu=0; - double s1u=0; - double sxu=0; - double syu=0; - double sxyu=0; - - opt_utilities::default_data_set<double,double> ds1; - for(size_t i=0;i<ds.size();++i) - { - double x=ds.get_data(i).get_x(); - double y=ds.get_data(i).get_y(); - double xl=ds.get_data(i).get_x_lower_err(); - double xu=ds.get_data(i).get_x_upper_err(); - double yl=ds.get_data(i).get_y_lower_err(); - double yu=ds.get_data(i).get_y_upper_err(); - - double new_x=shuffle_data(x, - xl, - xu); - double new_y=shuffle_data(y, - yl, - yu); - - ds1.add_data(data<double,double>(new_x,new_y, - yl/y*new_y, - yu/y*new_y, - xl/x*new_x, - xu/x*new_x)); - - x=new_x; - y=new_y; - if(x>Tb) - { - sxxl+=log(x)*log(x); - sxl+=log(x); - syl+=log(y); - sxyl+=log(y)*log(x); - s1l+=1; - } - else - { - sxxu+=log(x)*log(x); - sxu+=log(x); - syu+=log(y); - sxyu+=log(y)*log(x); - s1u+=1; - } - } - double Ml=sxxl*s1l-sxl*sxl; - double Mal=sxyl*s1l-syl*sxl; - double Mbl=sxxl*syl-sxl*sxyl; - double k0l=Mal/Ml; - double b0l=Mbl/Ml; - - double Mu=sxxu*s1u-sxu*sxu; - double Mau=sxyu*s1u-syu*sxu; - double Mbu=sxxu*syu-sxu*sxyu; - double k0u=Mau/Mu; - double b0u=Mbu/Mu; - - double gamma0l=k0l; - double gamma0u=k0u; - - double ampl0l=exp(b0l)*pow(Tb,gamma0l); - double ampl0u=exp(b0u)*pow(Tb,gamma0u);; - - fit.set_param_value("bpx",Tb); - fit.set_param_value("bpy",(ampl0l+ampl0u)/2); - fit.set_param_value("gamma1",gamma0l); - fit.set_param_value("gamma2",gamma0u); - - - fit.load_data(ds1); - - fit.fit(); - vector<double> p=fit.fit(); - for(size_t i=0;i<p.size();++i) - { - mean_p[i]+=p[i]; - mean_p2[i]+=p[i]*p[i]; - } - //cerr<<fit.get_param_value("gamma1")<<"\t"<<fit.get_param_value("gamma2")<<endl; - - } - vector<double> std_p(p.size()); - cerr<<endl; - for(size_t i=0;i<mean_p.size();++i) - { - mean_p[i]/=cnt; - mean_p2[i]/=cnt; - std_p[i]=std::sqrt(mean_p2[i]-mean_p[i]*mean_p[i]); - cout<<fit.get_param_info(i).get_name()<<"= "<<p[i]<<" +/- "<<std_p[i]<<endl; - } - std::cout<<"Num of sources:"<<ds.size()<<endl; -} diff --git a/mass_profile/fit_lt_pl.cpp b/mass_profile/fit_lt_pl.cpp deleted file mode 100644 index fa37249..0000000 --- a/mass_profile/fit_lt_pl.cpp +++ /dev/null @@ -1,239 +0,0 @@ -/* - Perform a double-beta density model fitting to the surface brightness data - Author: Junhua Gu - Last modified: 2011.01.01 - This code is distributed with no warrant -*/ - -//#define HAVE_X_ERROR -#include <iostream> -#include <sstream> -#include <fstream> -#include <models/pl1d.hpp> -#include <models/lin1d.hpp> -#include "statistics/chisq.hpp" -#include "statistics/leastsq.hpp" -#include <data_sets/default_data_set.hpp> -#include <methods/powell/powell_method.hpp> -#include <core/freeze_param.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; -const double pi=4*atan(1); -double std_norm_rand() -{ - double x=0; - double u=0; - double v=0; - - do - { - u=rand()/(double)RAND_MAX; - rand(); - v=rand()/(double)RAND_MAX; - - x=std::sqrt(-log(u))*cos(2*pi*v); - }while(isnan(x)); - return x; -} - -double shuffle_data(double xc,double xl,double xu) -{ - double result=0; - assert(!isnan(xc)); - assert(!isnan(xl)); - assert(!isnan(xu)); - if(std_norm_rand()>0) - { - result=xc-std::abs(std_norm_rand()*xl); - } - else - { - result=xc+std::abs(std_norm_rand()*xu); - } - assert(!isnan(result)); - return result; -} - -int main(int argc,char* argv[]) -{ - if(argc!=3) - { - cerr<<"Usage:"<<argv[0]<<" <a 5 column file with T -Terr +Terr L Lerr> <T lower limit>"<<endl; - return -1; - } - double T_lower_limit(atof(argv[2])); - ifstream ifs_data(argv[1]); - default_data_set<double,double> ds; - ofstream ofs_result("l-t_result.qdp"); - ofs_result<<"read terr 1"<<endl; - ofs_result<<"read serr 2"<<endl; - ofs_result<<"skip single"<<endl; - ofs_result<<"log"<<endl; - //ofs_result<<"li on 2"<<endl; - ofs_result<<"time off"<<endl; - ofs_result<<"la f"<<endl; - ofs_result<<"la x temperature (keV)"<<endl; - ofs_result<<"la y Luminosity (10\\u43\\d erg s\\u-1\\d)"<<endl; - double yunit=1E43; - double sxx=0; - double s1=0; - double sx=0; - double sy=0; - double sxy=0; - bool is_first_nonono=true; - for(;;) - { - double T,Tl,Tu; - double L,Lerr; - std::string line; - getline(ifs_data,line); - - - if(!ifs_data.good()) - { - break; - } - line+=" "; - istringstream iss(line); - - if(line[0]=='#') - { - if(!is_first_nonono) - { - ofs_result<<"no no no"<<endl; - } - else - { - is_first_nonono=false; - } - continue; - } - - iss>>T>>Tl>>Tu>>L>>Lerr; - //std::cerr<<L<<"\t"<<Lerr<<endl; - if(!iss.good()) - { - continue; - } - if(T<T_lower_limit||L<0) - { - continue; - } - if(Lerr<L*.1) - { - Lerr=L*.1; - } - double Ll=Lerr; - double Lu=Lerr; - Tl=std::abs(Tl); - Tu=std::abs(Tu); - Ll=std::abs(Ll); - Lu=std::abs(Lu); - ofs_result<<T<<"\t"<<-std::abs(Tl)<<"\t"<<+std::abs(Tu)<<"\t"<<L/yunit<<"\t"<<std::abs(Lerr)/yunit<<endl; - double x=log(T); - double y=log(L); - double xu=log(T+Tu)-log(T); - double xl=log(T-Tl)-log(T); - - double yu=log(L+Lu)-log(L); - double yl=log(L-Ll)-log(L); - if(isnan(x)||isnan(y)||isnan(yl)||isnan(yu)|| - isnan(xl)||isnan(xu)) - { - std::cerr<<"one data with error > data, skipped"<<endl; - std::cerr<<line<<endl; - continue; - } - sxx+=x*x; - sx+=x; - sy+=y; - sxy+=y*x; - s1+=1; - data<double,double> d(x,y,std::abs(yl),std::abs(yu), - std::abs(xl),std::abs(xu)); - ds.add_data(d); - } - - double M=sxx*s1-sx*sx; - double Ma=sxy*s1-sy*sx; - double Mb=sxx*sy-sx*sxy; - double k0=Ma/M; - double b0=Mb/M; - - ofs_result<<"no no no"<<endl; - fitter<double,double,vector<double>,double,std::string> fit; - fit.set_opt_method(powell_method<double,vector<double> >()); - fit.set_statistic(chisq<double,double,vector<double>,double,std::string>()); - //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>()); - fit.set_model(lin1d<double>()); - fit.load_data(ds); - - cerr<<"k0="<<k0<<endl; - cerr<<"b0="<<b0<<endl; - cerr<<"Ampl0="<<exp(b0)<<endl; - cerr<<"gamma0="<<k0<<endl; - fit.set_param_value("k",k0); - fit.set_param_value("b",b0); - fit.fit(); - std::vector<double> p=fit.fit(); - for(double i=.5;i<12;i*=1.01) - { - ofs_result<<i<<"\t0\t0\t"<<exp(fit.eval_model_raw(log(i),p))/yunit<<"\t0\n"; - } - - - double mean_A=0; - double mean_A2=0; - double mean_g=0; - double mean_g2=0; - int cnt=0; - for(int n=0;n<100;++n) - { - ++cnt; - cerr<<"."; - opt_utilities::default_data_set<double,double> ds1; - for(size_t i=0;i<ds.size();++i) - { - double new_x=shuffle_data(ds.get_data(i).get_x(), - ds.get_data(i).get_x_lower_err(), - ds.get_data(i).get_x_upper_err()); - double new_y=shuffle_data(ds.get_data(i).get_y(), - ds.get_data(i).get_y_lower_err(), - ds.get_data(i).get_y_upper_err()); - ds1.add_data(data<double,double>(new_x,new_y, - ds.get_data(i).get_y_lower_err(), - ds.get_data(i).get_y_upper_err(), - ds.get_data(i).get_y_lower_err(), - ds.get_data(i).get_y_upper_err())); - //cerr<<new_x<<"\t"<<new_y<<endl; - } - fit.load_data(ds1); - - fit.fit(); - double k=fit.get_param_value("k"); - double b=fit.get_param_value("b"); - double A=exp(b); - double g=k; - mean_A+=A; - mean_A2+=A*A; - mean_g+=g; - mean_g2+=g*g; - } - std::cerr<<endl; - mean_A/=cnt; - mean_A2/=cnt; - mean_g/=cnt; - mean_g2/=cnt; - double std_A=std::sqrt(mean_A2-mean_A*mean_A); - double std_g=std::sqrt(mean_g2-mean_g*mean_g); - - std::cerr<<"L=L0*T^gamma"<<endl; - std::cout<<"L0= "<<exp(p[1])<<"+/-"<<std_A<<endl; - std::cout<<"gamma= "<<p[0]<<"+/-"<<std_g<<endl; - std::cout<<"Num of sources:"<<ds.size()<<endl; - -} diff --git a/mass_profile/fit_mt_bpl.cpp b/mass_profile/fit_mt_bpl.cpp deleted file mode 100644 index e956f4f..0000000 --- a/mass_profile/fit_mt_bpl.cpp +++ /dev/null @@ -1,350 +0,0 @@ -/* - Perform a double-beta density model fitting to the surface brightness data - Author: Junhua Gu - Last modified: 2011.01.01 - This code is distributed with no warrant -*/ - -//#define HAVE_X_ERROR -#include <iomanip> -#include <iostream> -#include <sstream> -#include <fstream> -#include <models/bpl1d.hpp> -#include <models/lin1d.hpp> -#include "statistics/chisq.hpp" -#include "statistics/logchisq.hpp" -#include "statistics/leastsq.hpp" -#include <data_sets/default_data_set.hpp> -#include <methods/powell/powell_method.hpp> -#include <methods/gsl_simplex/gsl_simplex.hpp> -#include <core/freeze_param.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; -const double pi=4*atan(1); -double std_norm_rand() -{ - double u=0; - double v=0; - while(u<=0||v<=0) - { - u=rand()/(double)RAND_MAX; - rand(); - v=rand()/(double)RAND_MAX; - } - double x=std::sqrt(-log(u))*cos(2*pi*v); - return x; -} - -double shuffle_data(double xc,double xl,double xu) -{ - double lxc=log(xc); - double lxl=log(xc-xl)-log(xc); - double lxu=log(xc+xu)-log(xc); - - if(std_norm_rand()>0) - { - double result=std::exp(lxc-std::abs(std_norm_rand()*lxl)); - return result; - } - else - { - double result=std::exp(lxc+std::abs(std_norm_rand()*lxu)); - return result; - } -} - -int main(int argc,char* argv[]) -{ - srand(time(0)); - if(argc!=4) - { - cerr<<"Usage:"<<argv[0]<<" <a 6 column file with T -Terr +Terr M -Merr +Merr> <initial broken temperature> <T lower limit>"<<endl; - return -1; - } - double T_lower_limit(atof(argv[3])); - double Tb=atof(argv[2]); - assert(Tb>0); - ifstream ifs_data(argv[1]); - default_data_set<double,double> ds; - ofstream ofs_result("m-t_bpl_result.qdp"); - ofs_result<<"read terr 1 2"<<endl; - ofs_result<<"skip single"<<endl; - ofs_result<<"log"<<endl; - //ofs_result<<"li on 2"<<endl; - ofs_result<<"time off"<<endl; - ofs_result<<"la f"<<endl; - ofs_result<<"la x temperature (keV)"<<endl; - ofs_result<<"la y mass (M\\d\\(2281)\\u)"<<endl; - double sxxl=0; - double s1l=0; - double sxl=0; - double syl=0; - double sxyl=0; - - double sxxu=0; - double s1u=0; - double sxu=0; - double syu=0; - double sxyu=0; - - bool is_first_nonono=true; - - for(;;) - { - double T,Tl,Tu; - double M,Ml,Mu; - std::string line; - getline(ifs_data,line); - //ifs_data>>T>>Tl>>Tu>>M>>Ml>>Mu; - if(!ifs_data.good()) - { - break; - } - line+=" "; - istringstream iss(line); - - if(line[0]=='#') - { - if(!is_first_nonono) - { - ofs_result<<"no no no"<<endl; - } - else - { - is_first_nonono=false; - } - continue; - } - iss>>T>>Tl>>Tu>>M>>Ml>>Mu; - //std::cerr<<L<<"\t"<<Lerr<<endl; - if(!iss.good()) - { - continue; - } - - if(T<T_lower_limit||M<0) - { - continue; - } - - if(std::abs(Mu)<M*.1||std::abs(Ml)<M*.1) - { - cerr<<"mass error less than 10%, skipped"<<endl; - cerr<<line<<endl; - continue; - } -#if 1 - if(std::abs(Tu)<.1||std::abs(Tl)<.1) - { - cerr<<"T error less than 10%, skipped"<<endl; - cerr<<line<<endl; - continue; - } -#endif - - if(std::abs(Mu)+std::abs(Ml)<M*.1) - { - double k=M*.1/(std::abs(Mu)+std::abs(Ml)); - Mu*=k; - Ml*=k; - } - if(std::abs(Ml)>std::abs(M)) - { - continue; - } - Tl=std::abs(Tl); - Tu=std::abs(Tu); - Ml=std::abs(Ml); - Mu=std::abs(Mu); - ofs_result<<T<<"\t"<<-std::abs(Tl)<<"\t"<<+std::abs(Tu)<<"\t"<<M<<"\t"<<-std::abs(Ml)<<"\t"<<+std::abs(Mu)<<endl; - double x=(T); - double y=(M); - double xu=Tu; - double xl=Tl; - - double yu=Mu; - double yl=Ml; - if(T>Tb) - { - sxxl+=log(x)*log(x); - sxl+=log(x); - syl+=log(y); - sxyl+=log(y)*log(x); - s1l+=1; - } - else - { - sxxu+=log(x)*log(x); - sxu+=log(x); - syu+=log(y); - sxyu+=log(y)*log(x); - s1u+=1; - } - data<double,double> d(x,y,std::abs(yl),std::abs(yu), - std::abs(xl),std::abs(xu)); - ds.add_data(d); - } - - double Ml=sxxl*s1l-sxl*sxl; - double Mal=sxyl*s1l-syl*sxl; - double Mbl=sxxl*syl-sxl*sxyl; - double k0l=Mal/Ml; - double b0l=Mbl/Ml; - - double Mu=sxxu*s1u-sxu*sxu; - double Mau=sxyu*s1u-syu*sxu; - double Mbu=sxxu*syu-sxu*sxyu; - double k0u=Mau/Mu; - double b0u=Mbu/Mu; - - double gamma0l=k0l; - double gamma0u=k0u; - - double ampl0l=exp(b0l)*pow(Tb,gamma0l); - double ampl0u=exp(b0u)*pow(Tb,gamma0u);; - - - - ofs_result<<"no no no"<<endl; - fitter<double,double,vector<double>,double,std::string> fit; - fit.set_opt_method(powell_method<double,vector<double> >()); - - fit.set_statistic(logchisq<double,double,vector<double>,double,std::string>()); - //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>()); - fit.set_model(bpl1d<double>()); - fit.load_data(ds); - - cerr<<"k0l="<<k0l<<endl; - cerr<<"k0u="<<k0u<<endl; - cerr<<"Ampl0="<<(ampl0l+ampl0u)/2<<endl; - - fit.set_param_value("bpx",Tb); - fit.set_param_value("bpy",(ampl0l+ampl0u)/2); - fit.set_param_value("gamma1",gamma0l); - fit.set_param_value("gamma2",gamma0u); - - fit.fit(); - //fit.set_opt_method(gsl_simplex<double,vector<double> >()); - fit.fit(); - std::vector<double> p=fit.fit(); - Tb=fit.get_param_value("bpx"); - //std::cout<<"chi="<<fit.get_statistic().eval(p)<<std::endl; - for(double i=.5;i<12;i*=1.01) - { - ofs_result<<i<<"\t0\t0\t"<<fit.eval_model_raw(i,p)<<"\t0\t0\n"; - } - - - std::vector<double> mean_p(p.size()); - std::vector<double> mean_p2(p.size()); - int cnt=0; - for(int n=0;n<100;++n) - { - ++cnt; - cerr<<"."; - double sxxl=0; - double s1l=0; - double sxl=0; - double syl=0; - double sxyl=0; - - double sxxu=0; - double s1u=0; - double sxu=0; - double syu=0; - double sxyu=0; - - opt_utilities::default_data_set<double,double> ds1; - for(size_t i=0;i<ds.size();++i) - { - double x=ds.get_data(i).get_x(); - double y=ds.get_data(i).get_y(); - double xl=ds.get_data(i).get_x_lower_err(); - double xu=ds.get_data(i).get_x_upper_err(); - double yl=ds.get_data(i).get_y_lower_err(); - double yu=ds.get_data(i).get_y_upper_err(); - - double new_x=shuffle_data(x, - xl, - xu); - double new_y=shuffle_data(y, - yl, - yu); - - ds1.add_data(data<double,double>(new_x,new_y, - yl/y*new_y, - yu/y*new_y, - xl/x*new_x, - xu/x*new_x)); - - x=new_x; - y=new_y; - if(x>Tb) - { - sxxl+=log(x)*log(x); - sxl+=log(x); - syl+=log(y); - sxyl+=log(y)*log(x); - s1l+=1; - } - else - { - sxxu+=log(x)*log(x); - sxu+=log(x); - syu+=log(y); - sxyu+=log(y)*log(x); - s1u+=1; - } - } - double Ml=sxxl*s1l-sxl*sxl; - double Mal=sxyl*s1l-syl*sxl; - double Mbl=sxxl*syl-sxl*sxyl; - double k0l=Mal/Ml; - double b0l=Mbl/Ml; - - double Mu=sxxu*s1u-sxu*sxu; - double Mau=sxyu*s1u-syu*sxu; - double Mbu=sxxu*syu-sxu*sxyu; - double k0u=Mau/Mu; - double b0u=Mbu/Mu; - - double gamma0l=k0l; - double gamma0u=k0u; - - double ampl0l=exp(b0l)*pow(Tb,gamma0l); - double ampl0u=exp(b0u)*pow(Tb,gamma0u);; - - fit.set_param_value("bpx",Tb); - fit.set_param_value("bpy",(ampl0l+ampl0u)/2); - fit.set_param_value("gamma1",gamma0l); - fit.set_param_value("gamma2",gamma0u); - - - fit.load_data(ds1); - - fit.fit(); - vector<double> p=fit.fit(); - for(size_t i=0;i<p.size();++i) - { - mean_p[i]+=p[i]; - mean_p2[i]+=p[i]*p[i]; - } - //cerr<<fit.get_param_value("gamma1")<<"\t"<<fit.get_param_value("gamma2")<<endl; - - } - vector<double> std_p(p.size()); - cerr<<endl; - for(size_t i=0;i<mean_p.size();++i) - { - mean_p[i]/=cnt; - mean_p2[i]/=cnt; - std_p[i]=std::sqrt(mean_p2[i]-mean_p[i]*mean_p[i]); - cout<<fit.get_param_info(i).get_name()<<"= "<<p[i]<<" +/- "<<std_p[i]<<endl; - } - std::cout<<"Num of sources:"<<ds.size()<<endl; -} diff --git a/mass_profile/fit_mt_pl.cpp b/mass_profile/fit_mt_pl.cpp deleted file mode 100644 index bfccaf4..0000000 --- a/mass_profile/fit_mt_pl.cpp +++ /dev/null @@ -1,265 +0,0 @@ -/* - Perform a double-beta density model fitting to the surface brightness data - Author: Junhua Gu - Last modified: 2011.01.01 - This code is distributed with no warrant -*/ - -//#define HAVE_X_ERROR -#include <iomanip> -#include <iostream> -#include <sstream> -#include <fstream> -#include <models/pl1d.hpp> -#include <models/lin1d.hpp> -#include "statistics/chisq.hpp" -#include "statistics/leastsq.hpp" -#include "statistics/robust_chisq.hpp" -#include <data_sets/default_data_set.hpp> -#include <methods/powell/powell_method.hpp> -#include <core/freeze_param.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; -const double pi=4*atan(1); -double std_norm_rand() -{ - double u=0; - double v=0; - while(u<=0||v<=0) - { - u=rand()/(double)RAND_MAX; - rand(); - v=rand()/(double)RAND_MAX; - } - double x=std::sqrt(-log(u))*cos(2*pi*v); - return x; -} - -double shuffle_data(double xc,double xl,double xu) -{ - if(std_norm_rand()>0) - { - double result=xc-std::abs(std_norm_rand()*xl); - return result; - } - else - { - double result= xc+std::abs(std_norm_rand()*xu); - return result; - } -} - -int main(int argc,char* argv[]) -{ - if(argc!=3) - { - cerr<<"Usage:"<<argv[0]<<" <a 6 column file with T -Terr +Terr M -Merr +Merr> <lower T limit>"<<endl; - return -1; - } - double T_lower_limit(atof(argv[2])); - ifstream ifs_data(argv[1]); - default_data_set<double,double> ds; - ofstream ofs_result("m-t_result.qdp"); - ofs_result<<"read terr 1 2"<<endl; - ofs_result<<"skip single"<<endl; - ofs_result<<"log"<<endl; - //ofs_result<<"li on 2"<<endl; - ofs_result<<"time off"<<endl; - ofs_result<<"la f"<<endl; - ofs_result<<"la x temperature (keV)"<<endl; - ofs_result<<"la y mass (M\\dsun\\u)"<<endl; - double sxx=0; - double s1=0; - double sx=0; - double sy=0; - double sxy=0; - bool is_first_nonono=true; - - for(;;) - { - double T,Tl,Tu; - double M,Ml,Mu; - std::string line; - getline(ifs_data,line); - //ifs_data>>T>>Tl>>Tu>>M>>Ml>>Mu; - if(!ifs_data.good()) - { - break; - } - line+=" "; - istringstream iss(line); - - if(line[0]=='#') - { - if(!is_first_nonono) - { - ofs_result<<"no no no"<<endl; - } - else - { - is_first_nonono=false; - } - continue; - } - iss>>T>>Tl>>Tu>>M>>Ml>>Mu; - //std::cerr<<L<<"\t"<<Lerr<<endl; - if(!iss.good()) - { - continue; - } - - if(T<T_lower_limit||M<0) - { - continue; - } - if(std::abs(Mu)<M*.1||std::abs(Ml)<M*.1) - { - cerr<<"mass error less than 10%, skipped"<<endl; - cerr<<line<<endl; - continue; - } -#if 1 - if(std::abs(Tu)<.1||std::abs(Tl)<.1) - { - cerr<<"T error less than 10%, skipped"<<endl; - cerr<<line<<endl; - continue; - } -#endif - if(std::abs(Mu)+std::abs(Ml)<M*.1) - { - double k=M*.1/(std::abs(Mu)+std::abs(Ml)); - Mu*=k; - Ml*=k; - } - Tl=std::abs(Tl); - Tu=std::abs(Tu); - Ml=std::abs(Ml); - Mu=std::abs(Mu); - ofs_result<<T<<"\t"<<-std::abs(Tl)<<"\t"<<+std::abs(Tu)<<"\t"<<M<<"\t"<<-std::abs(Ml)<<"\t"<<+std::abs(Mu)<<endl; - double x=log(T); - double y=log(M); - double xu=log(T+Tu)-log(T); - double xl=log(T-Tl)-log(T); - - double yu=log(M+Mu)-log(M); - double yl=log(M-Ml)-log(M); - if(isnan(x)||isnan(y)||isnan(yl)||isnan(yu)|| - isnan(xl)||isnan(xu)) - { - std::cerr<<"one data with error > data, skipped"<<endl; - std::cerr<<line<<endl; - continue; - } - sxx+=x*x; - sx+=x; - sy+=y; - sxy+=y*x; - s1+=1; - data<double,double> d(x,y,std::abs(yl),std::abs(yu), - std::abs(xl),std::abs(xu)); - ds.add_data(d); - } - - double M=sxx*s1-sx*sx; - double Ma=sxy*s1-sy*sx; - double Mb=sxx*sy-sx*sxy; - double k0=Ma/M; - double b0=Mb/M; - - ofs_result<<"no no no"<<endl; - fitter<double,double,vector<double>,double,std::string> fit; - fit.set_opt_method(powell_method<double,vector<double> >()); - fit.set_statistic(chisq<double,double,vector<double>,double,std::string>()); - //fit.set_statistic(robust_chisq<double,double,vector<double>,double,std::string>()); - //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>()); - fit.set_model(lin1d<double>()); - fit.load_data(ds); - - cerr<<"k0="<<k0<<endl; - cerr<<"b0="<<b0<<endl; - cerr<<"Ampl0="<<exp(b0)<<endl; - cerr<<"gamma0="<<k0<<endl; - fit.set_param_value("k",k0); - fit.set_param_value("b",b0); - std::vector<double> p=fit.get_all_params(); - std::cout<<"chi="<<fit.get_statistic().eval(p)<<std::endl; - fit.fit(); - fit.fit(); - p=fit.fit(); - - std::cout<<"chi="<<fit.get_statistic().eval(p)<<std::endl; - for(double i=.5;i<12;i*=1.01) - { - ofs_result<<i<<"\t0\t0\t"<<exp(fit.eval_model_raw(log(i),p))<<"\t0\t0\n"; - } - - ofstream ofs_resid("resid.qdp"); - ofs_resid<<"read terr 1 2 3"<<endl; - ofs_resid<<"skip single"<<endl; - ofs_resid<<"ma 3 on 1"<<endl; - ofs_resid<<"log x"<<endl; - for(size_t i=0;i<ds.size();++i) - { - double x=ds.get_data(i).get_x(); - double y=ds.get_data(i).get_y(); - //double xe1=-ds.get_data(i).get_x_lower_err()*0; - //double xe2=ds.get_data(i).get_x_upper_err()*0; - double ye1=-ds.get_data(i).get_y_lower_err(); - double ye2=ds.get_data(i).get_y_upper_err(); - ofs_resid<<exp(x)<<"\t"<<0<<"\t"<<0<<"\t"<<y-fit.eval_model_raw(x,p)<<"\t"<<ye1<<"\t"<<ye2<<"\t"<<"0\t0\t0"<<endl; - } - double mean_A=0; - double mean_A2=0; - double mean_g=0; - double mean_g2=0; - int cnt=0; - for(int n=0;n<100;++n) - { - ++cnt; - cerr<<"."; - opt_utilities::default_data_set<double,double> ds1; - for(size_t i=0;i<ds.size();++i) - { - double new_x=shuffle_data(ds.get_data(i).get_x(), - ds.get_data(i).get_x_lower_err(), - ds.get_data(i).get_x_upper_err()); - double new_y=shuffle_data(ds.get_data(i).get_y(), - ds.get_data(i).get_y_lower_err(), - ds.get_data(i).get_y_upper_err()); - ds1.add_data(data<double,double>(new_x,new_y, - ds.get_data(i).get_y_lower_err(), - ds.get_data(i).get_y_upper_err(), - ds.get_data(i).get_y_lower_err(), - ds.get_data(i).get_y_upper_err())); - } - fit.load_data(ds1); - - fit.fit(); - double k=fit.get_param_value("k"); - double b=fit.get_param_value("b"); - double A=exp(b); - double g=k; - mean_A+=A; - mean_A2+=A*A; - mean_g+=g; - mean_g2+=g*g; - } - std::cerr<<endl; - mean_A/=cnt; - mean_A2/=cnt; - mean_g/=cnt; - mean_g2/=cnt; - double std_A=std::sqrt(mean_A2-mean_A*mean_A); - double std_g=std::sqrt(mean_g2-mean_g*mean_g); - - std::cerr<<"M=M0*T^gamma"<<endl; - std::cout<<"M0= "<<exp(p[1])<<"+/-"<<std_A<<endl; - std::cout<<"gamma= "<<p[0]<<"+/-"<<std_g<<endl; - std::cout<<"Num of sources:"<<ds.size()<<endl; - -} |