aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--mass_profile/beta.hpp4
-rw-r--r--mass_profile/calc_lx.cpp16
-rw-r--r--mass_profile/calc_lx_beta.cpp79
-rw-r--r--mass_profile/calc_lx_dbeta.cpp103
-rw-r--r--mass_profile/chisq.hpp26
-rw-r--r--mass_profile/dbeta.hpp10
-rw-r--r--mass_profile/dump_fit_qdp.cpp14
-rw-r--r--mass_profile/fit_beta_sbp.cpp40
-rw-r--r--mass_profile/fit_dbeta_sbp.cpp55
-rw-r--r--mass_profile/fit_direct_beta.cpp10
-rw-r--r--mass_profile/fit_lt_bpl.cpp42
-rw-r--r--mass_profile/fit_lt_pl.cpp22
-rw-r--r--mass_profile/fit_mt_bpl.cpp44
-rw-r--r--mass_profile/fit_mt_pl.cpp18
-rw-r--r--mass_profile/fit_nfw_mass.cpp22
-rw-r--r--mass_profile/fit_nfw_sbp.cpp34
-rw-r--r--mass_profile/fit_wang2012_model.cpp20
-rw-r--r--mass_profile/nfw_ne.hpp20
-rw-r--r--mass_profile/projector.hpp8
-rw-r--r--mass_profile/vchisq.hpp32
20 files changed, 300 insertions, 319 deletions
diff --git a/mass_profile/beta.hpp b/mass_profile/beta.hpp
index 4e6264c..6ae70ac 100644
--- a/mass_profile/beta.hpp
+++ b/mass_profile/beta.hpp
@@ -13,7 +13,7 @@ namespace opt_utilities
{
this->push_param_info(param_info<std::vector<T>,std::string>("n0",1,0,1E99));
this->push_param_info(param_info<std::vector<T>,std::string>("beta",.66,0,1E99));
- this->push_param_info(param_info<std::vector<T>,std::string>("rc",100,0,1E99));
+ this->push_param_info(param_info<std::vector<T>,std::string>("rc",100,0,1E99));
}
public:
@@ -30,7 +30,7 @@ namespace opt_utilities
T rc=p[2];
std::vector<T> result(x.size()-1);
- for(int i=1;i<x.size();++i)
+ for(size_t i=1;i<x.size();++i)
{
T xi=(x[i]+x[i-1])/2;
T yi=0;
diff --git a/mass_profile/calc_lx.cpp b/mass_profile/calc_lx.cpp
index e5cff62..b1137ff 100644
--- a/mass_profile/calc_lx.cpp
+++ b/mass_profile/calc_lx.cpp
@@ -61,7 +61,7 @@ int main(int argc,char* argv[])
//perform a 1-beta fitting////
//////////////////////////////
fitter<double,double,vector<double>,double,string> f;
-
+
f.set_statistic(chisq<double,double,vector<double>,double,string>());
f.set_opt_method(powell_method<double,vector<double> >());
f.set_model(beta1d<double>());
@@ -74,7 +74,7 @@ int main(int argc,char* argv[])
double rmax=f.get_data_set().get_data(f.get_data_set().size()-1).get_x();
ofstream lx_fit_result("lx_fit_result.qdp");
lx_fit_result<<"read terr 1 2\nskip single\n";
- for(int i=0;i<f.get_data_set().size();++i)
+ for(size_t i=0;i<f.get_data_set().size();++i)
{
lx_fit_result<<f.get_data_set().get_data(i).get_x()<<"\t"<<
-abs(f.get_data_set().get_data(i).get_x_lower_err())<<"\t"<<
@@ -89,13 +89,13 @@ int main(int argc,char* argv[])
{
lx_fit_result<<i<<"\t0\t0\t"<<f.eval_model(i,f.get_all_params())<<"\t0\t0"<<endl;
}
-
- for(int i=0;i<f.get_num_params();++i)
+
+ for(size_t i=0;i<f.get_num_params();++i)
{
cerr<<f.get_param_info(i).get_name()<<"="<<
f.get_param_info(i).get_value()<<endl;
}
-
+
vector<double> p=f.get_all_params();
double r500kpc=atof(argv[4]);
@@ -131,7 +131,7 @@ int main(int argc,char* argv[])
sum_flux+=cnt*ratio;
}
double lx=sum_flux*4*pi*d_l_cm*d_l_cm;
-
+
double l_mean=0;
double l2_mean=0;
double cnt=0;
@@ -140,7 +140,7 @@ int main(int argc,char* argv[])
cerr<<".";
opt_utilities::default_data_set<double,double> ds(dl.get_data_set());
opt_utilities::default_data_set<double,double> ds1;
- for(int i=0;i<ds.size();++i)
+ for(size_t i=0;i<ds.size();++i)
{
double yc=ds.get_data(i).get_y();
double ye=(std::abs(ds.get_data(i).get_y_lower_err())+std::abs(ds.get_data(i).get_y_lower_err()))/2;
@@ -159,7 +159,7 @@ int main(int argc,char* argv[])
//cout<<f.get_param_value("beta")<<endl;
double sum_cnt=0;
double sum_flux=0;
-
+
for(double r=0;r<r500pixel;r+=dr)
{
double r1=r<.2*r500pixel?.2*r500pixel:r;
diff --git a/mass_profile/calc_lx_beta.cpp b/mass_profile/calc_lx_beta.cpp
index 5bb2570..aae8b2b 100644
--- a/mass_profile/calc_lx_beta.cpp
+++ b/mass_profile/calc_lx_beta.cpp
@@ -30,20 +30,8 @@ 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;
+ return abs(n0)*pow(1+r*r/rc/rc,-3./2.*abs(beta));
}
@@ -61,7 +49,7 @@ public:
{
return spl.get_value(x);
}
-
+
//we need this function, when this object is performing a clone of itself
spline_func_obj* do_clone()const
{
@@ -74,7 +62,7 @@ public:
{
spl.push_point(x,y);
}
-
+
//before getting the intepolated value, the spline should be initialzied by calling this function
void gen_spline()
{
@@ -93,9 +81,9 @@ int main(int argc,char* argv[])
ifstream cfg_file(argv[1]);
assert(cfg_file.is_open());
cfg_map cfg=parse_cfg_file(cfg_file);
-
+
//check the existence of following parameters
-
+
const double z=cfg.z;
//initialize the radius list, sbp list and sbp error list
@@ -151,7 +139,7 @@ int main(int argc,char* argv[])
{
rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
}
-
+
cerr<<"rmin="<<rmin<<endl;
std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
radii_tmp.resize(radii.size());
@@ -184,7 +172,7 @@ int main(int argc,char* argv[])
for(ifstream ifs(cfg.cfunc_file.c_str());;)
{
assert(ifs.is_open());
- double x,y,y1,y2;
+ double x,y;
ifs>>x>>y;
if(!ifs.good())
{
@@ -199,7 +187,7 @@ int main(int argc,char* argv[])
cf.add_point(x,y);//change with source
}
cf.gen_spline();
-
+
//read temperature profile data
spline_func_obj Tprof;
int tcnt=0;
@@ -224,7 +212,7 @@ int main(int argc,char* argv[])
Tprof.gen_spline();
-
+
default_data_set<std::vector<double>,std::vector<double> > ds;
ds.add_data(data<std::vector<double>,std::vector<double> >(radii,sbps,sbpe,sbpe,radii,radii));
@@ -247,11 +235,12 @@ int main(int argc,char* argv[])
//optimization method
f.set_opt_method(powell_method<double,std::vector<double> >());
//initialize the initial values
+ /*
double n0=0;
- //double beta=atof(arg_map["beta"].c_str());
double beta=0;
double rc=0;
- double bkg=0;
+ */
+ double bkg_level=0;
for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin();
i!=cfg.param_map.end();++i)
@@ -280,12 +269,14 @@ int main(int argc,char* argv[])
f.fit();
f.fit();
std::vector<double> 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(int i=0;i<f.get_num_params();++i)
+ for(size_t i=0;i<f.get_num_params();++i)
{
if(f.get_param_info(i).get_name()=="rc")
{
@@ -322,7 +313,7 @@ int main(int argc,char* argv[])
ofs_sbp<<"color 3 on 3"<<endl;
ofs_sbp<<"color 4 on 4"<<endl;
ofs_sbp<<"color 5 on 5"<<endl;
-
+
ofs_sbp<<"win 1"<<endl;
ofs_sbp<<"yplot 1 2 3 4 5"<<endl;
ofs_sbp<<"loc 0 0 1 1"<<endl;
@@ -356,7 +347,7 @@ int main(int argc,char* argv[])
ofs_sbp<<"color 4 on 3"<<endl;
ofs_sbp<<"color 5 on 4"<<endl;
//ofs_sbp<<"ma 1 on 2"<<endl;
-
+
ofs_sbp<<"win 1"<<endl;
ofs_sbp<<"yplot 1 2 3 4"<<endl;
ofs_sbp<<"loc 0 0 1 1"<<endl;
@@ -377,7 +368,7 @@ int main(int argc,char* argv[])
}
// cout<<sbps_all.size()<<"\t"<<sbps.size()<<"\t"<<sbps_inner_cut_size<<endl;
- for(int i=1;i<sbps_all.size();++i)
+ for(size_t i=1;i<sbps_all.size();++i)
{
double x=(radii_all[i]+radii_all[i-1])/2;
double y=sbps_all[i-1];
@@ -395,7 +386,7 @@ int main(int argc,char* argv[])
}
}
ofs_sbp<<"no no no"<<endl;
- for(int i=sbps_inner_cut_size;i<sbps_all.size();++i)
+ for(size_t i=sbps_inner_cut_size;i<sbps_all.size();++i)
{
double x=(radii_all[i]+radii_all[i-1])/2;
double ym=mv[i-1];
@@ -403,8 +394,8 @@ int main(int argc,char* argv[])
}
ofs_sbp<<"no no no"<<endl;
//bkg level
- double bkg_level=abs(f.get_param_value("bkg"));
- for(int i=0;i<sbps_all.size();++i)
+ bkg_level=abs(f.get_param_value("bkg"));
+ for(size_t i=0;i<sbps_all.size();++i)
{
double x=(radii_all[i]+radii_all[i-1])/2;
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
@@ -420,23 +411,23 @@ int main(int argc,char* argv[])
}
//resid
ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
+ //double y=sbps[i-1];
+ //double ye=sbpe[i-1];
double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
}
//zero level of resid
ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
- double ym=mv[i-1];
+ //double y=sbps[i-1];
+ //double ye=sbpe[i-1];
+ //double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl;
}
@@ -457,13 +448,13 @@ int main(int argc,char* argv[])
p.back()=0;
radii.clear();
double rout=atof(argv[2])*kpc;
-
+
for(double r=0;r<rout;r+=1*kpc)//step size=1kpc
{
double r_pix=r/cm_per_pixel;
radii.push_back(r_pix);
}
-
+
double Da=cm_per_pixel/(.492/3600./180.*pi);
double Dl=Da*(1+z)*(1+z);
cout<<"dl="<<Dl/kpc<<endl;
@@ -481,19 +472,19 @@ int main(int argc,char* argv[])
break;
}
//cerr<<x<<"\t"<<y<<endl;
-
+
cf_bolo_erg.add_point(x,y);//change with source
}
cf_bolo_erg.gen_spline();
projector<double>& pj=dynamic_cast<projector<double>&>(f.get_model());
pj.attach_cfunc(cf_bolo_erg);
-
-
-
+
+
+
mv=f.eval_model_raw(radii,p);
double flux_erg=0;
- for(int i=0;i<radii.size()-1;++i)
+ for(size_t i=0;i<radii.size()-1;++i)
{
double S=pi*(radii[i+1]+radii[i])*(radii[i+1]-radii[i]);
flux_erg+=S*mv[i];
diff --git a/mass_profile/calc_lx_dbeta.cpp b/mass_profile/calc_lx_dbeta.cpp
index f105869..f8f2c2a 100644
--- a/mass_profile/calc_lx_dbeta.cpp
+++ b/mass_profile/calc_lx_dbeta.cpp
@@ -26,20 +26,8 @@ double dbeta_func(double r,
double n01,double rc1,double beta1,
double n02,double rc2,double beta2)
{
-
- return abs(n01)*pow(1+r*r/rc1/rc1,-3./2.*abs(beta1))+abs(n02)*pow(1+r*r/rc2/rc2,-3./2.*abs(beta2));
-}
-
- //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;
+ return abs(n01)*pow(1+r*r/rc1/rc1,-3./2.*abs(beta1))+abs(n02)*pow(1+r*r/rc2/rc2,-3./2.*abs(beta2));
}
@@ -57,7 +45,7 @@ public:
{
return spl.get_value(x);
}
-
+
//we need this function, when this object is performing a clone of itself
spline_func_obj* do_clone()const
{
@@ -70,7 +58,7 @@ public:
{
spl.push_point(x,y);
}
-
+
//before getting the intepolated value, the spline should be initialzied by calling this function
void gen_spline()
{
@@ -89,9 +77,9 @@ int main(int argc,char* argv[])
ifstream cfg_file(argv[1]);
assert(cfg_file.is_open());
cfg_map cfg=parse_cfg_file(cfg_file);
-
+
//check the existence of following parameters
-
+
const double z=cfg.z;
//initialize the radius list, sbp list and sbp error list
@@ -147,7 +135,7 @@ int main(int argc,char* argv[])
{
rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
}
-
+
cerr<<"rmin="<<rmin<<endl;
std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
radii_tmp.resize(radii.size());
@@ -180,7 +168,7 @@ int main(int argc,char* argv[])
for(ifstream ifs(cfg.cfunc_file.c_str());;)
{
assert(ifs.is_open());
- double x,y,y1,y2;
+ double x,y;
ifs>>x>>y;
if(!ifs.good())
{
@@ -195,7 +183,7 @@ int main(int argc,char* argv[])
cf.add_point(x,y);//change with source
}
cf.gen_spline();
-
+
//read temperature profile data
spline_func_obj Tprof;
int tcnt=0;
@@ -220,10 +208,10 @@ int main(int argc,char* argv[])
Tprof.gen_spline();
-
+
default_data_set<std::vector<double>,std::vector<double> > ds;
ds.add_data(data<std::vector<double>,std::vector<double> >(radii,sbps,sbpe,sbpe,radii,radii));
-
+
//initial fitter
fitter<vector<double>,vector<double>,vector<double>,double> f;
f.load_data(ds);
@@ -255,7 +243,7 @@ int main(int argc,char* argv[])
//attach the cooling function
a.attach_cfunc(cf);
a.set_cm_per_pixel(cm_per_pixel);
-
+
f.set_model(a);
//chi^2 statistic
vchisq<double> c;
@@ -265,12 +253,15 @@ int main(int argc,char* argv[])
//optimization method
f.set_opt_method(powell_method<double,std::vector<double> >());
//initialize the initial values
- double n01=0;
+ /*
+ double =0;
double rc1=0;
double n02=0;
double rc2=0;
double beta=0;
- double bkg=0;
+ */
+ double bkg_level=0;
+
if(tie_beta)
{
f.set_param_value("beta",.7);
@@ -311,7 +302,7 @@ int main(int argc,char* argv[])
}
-
+
//perform the fitting, first freeze beta1, beta2, rc1, and rc2
if(tie_beta)
{
@@ -328,9 +319,9 @@ int main(int argc,char* argv[])
freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc2")
);
}
-
+
f.fit();
-
+
f.clear_param_modifier();
//then perform the fitting, freeze beta1 and beta2
@@ -341,9 +332,11 @@ int main(int argc,char* argv[])
//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");
@@ -359,10 +352,12 @@ int main(int argc,char* argv[])
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(int i=0;i<f.get_num_params();++i)
+ for(size_t i=0;i<f.get_num_params();++i)
{
if(f.get_param_info(i).get_name()=="rc1")
{
@@ -379,19 +374,19 @@ int main(int argc,char* argv[])
}
cerr<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
param_output<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-
+
//c.verbose(false);
//f.set_statistic(c);
//f.fit();
std::vector<double> p=f.get_all_params();
f.clear_param_modifier();
std::vector<double> mv=f.eval_model(radii,p);
-
-
+
+
ofstream ofs_sbp("lx_sbp_fit.qdp");
ofs_sbp<<"read serr 2"<<endl;
ofs_sbp<<"skip single"<<endl;
-
+
ofs_sbp<<"line on 2"<<endl;
ofs_sbp<<"line on 3"<<endl;
ofs_sbp<<"line on 4"<<endl;
@@ -401,7 +396,7 @@ int main(int argc,char* argv[])
ofs_sbp<<"ls 2 on 3"<<endl;
ofs_sbp<<"ls 2 on 4"<<endl;
ofs_sbp<<"ls 2 on 5"<<endl;
-
+
ofs_sbp<<"!LAB POS Y 4.00"<<endl;
@@ -423,27 +418,27 @@ int main(int argc,char* argv[])
ofs_sbp<<"log x"<<endl;
ofs_sbp<<"log y off"<<endl;
ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
double y=sbps[i-1];
double ye=sbpe[i-1];
- double ym=mv[i-1];
+ //double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<endl;
}
ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
+ //double y=sbps[i-1];
+ //double ye=sbpe[i-1];
double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<0<<endl;
}
//bkg
ofs_sbp<<"no no no"<<endl;
- double bkg_level=abs(f.get_param_value("bkg"));
- for(int i=0;i<sbps.size();++i)
+ bkg_level=abs(f.get_param_value("bkg"));
+ for(size_t i=0;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
@@ -466,22 +461,22 @@ int main(int argc,char* argv[])
}
//resid
ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
+ //double y=sbps[i-1];
+ //double ye=sbpe[i-1];
double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
}
//zero level in resid map
ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
- double ym=mv[i-1];
+ //double y=sbps[i-1];
+ //double ye=sbpe[i-1];
+ //double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl;
}
@@ -502,7 +497,7 @@ int main(int argc,char* argv[])
*/
p.back()=0;
-
+
radii.clear();
double rout=atof(argv[2])*kpc;
for(double r=0;r<rout;r+=1*kpc)//step size=1kpc
@@ -527,17 +522,17 @@ int main(int argc,char* argv[])
break;
}
//cerr<<x<<"\t"<<y<<endl;
-
+
cf_bolo_erg.add_point(x,y);//change with source
}
cf_bolo_erg.gen_spline();
-
+
projector<double>& pj=dynamic_cast<projector<double>&>(f.get_model());
pj.attach_cfunc(cf_bolo_erg);
-
+
mv=f.eval_model_raw(radii,p);
double flux_erg=0;
- for(int i=0;i<radii.size()-1;++i)
+ for(size_t i=0;i<radii.size()-1;++i)
{
double S=pi*(radii[i+1]+radii[i])*(radii[i+1]-radii[i]);
flux_erg+=S*mv[i];
diff --git a/mass_profile/chisq.hpp b/mass_profile/chisq.hpp
index 575f1c0..bdfadcb 100644
--- a/mass_profile/chisq.hpp
+++ b/mass_profile/chisq.hpp
@@ -46,7 +46,7 @@ namespace opt_utilities
bool verb;
bool limit_bound;
int n;
-
+
statistic<Ty,Tx,Tp,Ts,Tstr>* do_clone()const
{
// return const_cast<statistic<Ty,Tx,Tp>*>(this);
@@ -76,15 +76,15 @@ namespace opt_utilities
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(int i=0;i<p1.size();++i)
+ for(size_t i=0;i<p1.size();++i)
{
if(p1[i]>this->get_fitter().get_param_info(i).get_upper_limit()||
p1[i]<this->get_fitter().get_param_info(i).get_lower_limit())
@@ -111,7 +111,7 @@ namespace opt_utilities
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)
@@ -132,7 +132,7 @@ namespace opt_utilities
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(errx1<errx2)
{
@@ -170,7 +170,7 @@ namespace opt_utilities
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();
@@ -178,7 +178,7 @@ namespace opt_utilities
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);
@@ -192,7 +192,7 @@ namespace opt_utilities
if(n%display_interval==0)
{
cerr<<result<<"\t";
- for(int i=0;i<(int)get_size(p);++i)
+ for(size_t i=0;i<get_size(p);++i)
{
cerr<<get_element(p,i)<<",";
}
@@ -205,15 +205,13 @@ namespace opt_utilities
}
}
-
+
return result;
}
};
-
-
+
+
}
#endif
//EOF
-
-
diff --git a/mass_profile/dbeta.hpp b/mass_profile/dbeta.hpp
index 8c3a7c5..7246b44 100644
--- a/mass_profile/dbeta.hpp
+++ b/mass_profile/dbeta.hpp
@@ -24,7 +24,7 @@ namespace opt_utilities
this->push_param_info(param_info<std::vector<T>,std::string>("n02",1));
this->push_param_info(param_info<std::vector<T>,std::string>("beta2",.67));
this->push_param_info(param_info<std::vector<T>,std::string>("rc2",110));
-
+
}
public:
@@ -44,10 +44,10 @@ namespace opt_utilities
T beta2=p[4];
T rc2=p[5];
-
+
std::vector<T> result(x.size()-1);
- for(int i=1;i<x.size();++i)
+ for(size_t i=1;i<x.size();++i)
{
T xi=(x[i]+x[i-1])/2;
T yi=0;
@@ -70,7 +70,7 @@ namespace opt_utilities
this->push_param_info(param_info<std::vector<T>,std::string>("n02",1));
this->push_param_info(param_info<std::vector<T>,std::string>("rc2",110));
this->push_param_info(param_info<std::vector<T>,std::string>("beta",.67));
-
+
}
public:
@@ -92,7 +92,7 @@ namespace opt_utilities
T beta2=beta;
std::vector<T> result(x.size()-1);
- for(int i=1;i<x.size();++i)
+ for(size_t i=1;i<x.size();++i)
{
T xi=(x[i]+x[i-1])/2;
T yi=0;
diff --git a/mass_profile/dump_fit_qdp.cpp b/mass_profile/dump_fit_qdp.cpp
index 556edea..85307d1 100644
--- a/mass_profile/dump_fit_qdp.cpp
+++ b/mass_profile/dump_fit_qdp.cpp
@@ -10,13 +10,13 @@ namespace opt_utilities
os<<"la x \"radius (kpc)\""<<std::endl;
os<<"la y \"surface brightness (cts s\\u-1\\d pixel\\u-2\\d)\""<<std::endl;
os<<"li on 2"<<std::endl;
- for(int i=1;i<r.size();++i)
+ for(size_t i=1;i<r.size();++i)
{
os<<(r[i]+r[i-1])/2*cm_per_pixel/kpc<<"\t"<<(r[i]-r[i-1])/2*cm_per_pixel/kpc<<"\t"<<y[i-1]<<"\t"<<ye[i-1]<<std::endl;
}
os<<"no no no"<<std::endl;
std::vector<double> p=f.get_all_params();
- for(int i=1;i<r.size();++i)
+ for(size_t i=1;i<r.size();++i)
{
double x=(r[i]+r[i-1])/2;
os<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<f.eval_model_raw(x,p)<<"\t"<<0<<std::endl;
@@ -30,26 +30,26 @@ namespace opt_utilities
os<<"la x \"radius (kpc)\""<<std::endl;
os<<"la y \"density (cm\\u-3\\d)\""<<std::endl;
os<<"li on 2"<<std::endl;
-
- for(int i=1;i<r.size();++i)
+
+ for(size_t i=1;i<r.size();++i)
{
double x=(r[i]+r[i-1])/2;
double y=sbps[i-1];
double ye=sbpe[i-1];
os<<x*cm_per_pixel/kpc<<"\t0\t"<<y<<"\t"<<ye<<std::endl;
}
-
+
os<<"no no no"<<std::endl;
std::vector<double> p=f.get_all_params();
std::vector<double> mv=f.eval_model_raw(r,p);
- for(int i=1;i<r.size();++i)
+ for(size_t i=1;i<r.size();++i)
{
double x=(r[i]+r[i-1])/2;
double y=mv[i-1];
double ye=0;
os<<x*cm_per_pixel/kpc<<"\t0\t"<<y<<"\t"<<ye<<std::endl;
}
-
+
}
};
diff --git a/mass_profile/fit_beta_sbp.cpp b/mass_profile/fit_beta_sbp.cpp
index 41d6192..9db48cc 100644
--- a/mass_profile/fit_beta_sbp.cpp
+++ b/mass_profile/fit_beta_sbp.cpp
@@ -184,7 +184,7 @@ int main(int argc,char* argv[])
for(ifstream ifs(cfg.cfunc_file.c_str());;)
{
assert(ifs.is_open());
- double x,y,y1,y2;
+ double x,y;
ifs>>x>>y;
if(!ifs.good())
{
@@ -251,7 +251,7 @@ int main(int argc,char* argv[])
//double beta=atof(arg_map["beta"].c_str());
double beta=0;
double rc=0;
- double bkg=0;
+ double bkg_level=0;
for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin();
i!=cfg.param_map.end();++i)
@@ -285,7 +285,7 @@ int main(int argc,char* argv[])
beta=f.get_param_value("beta");
//output the datasets and fitting results
ofstream param_output("beta_param.txt");
- for(int i=0;i<f.get_num_params();++i)
+ for(size_t i=0;i<f.get_num_params();++i)
{
if(f.get_param_info(i).get_name()=="rc")
{
@@ -377,7 +377,7 @@ int main(int argc,char* argv[])
}
// cout<<sbps_all.size()<<"\t"<<sbps.size()<<"\t"<<sbps_inner_cut_size<<endl;
- for(int i=1;i<sbps_all.size();++i)
+ for(size_t i=1;i<sbps_all.size();++i)
{
double x=(radii_all[i]+radii_all[i-1])/2;
double y=sbps_all[i-1];
@@ -395,7 +395,7 @@ int main(int argc,char* argv[])
}
}
ofs_sbp<<"no no no"<<endl;
- for(int i=sbps_inner_cut_size;i<sbps_all.size();++i)
+ for(size_t i=sbps_inner_cut_size;i<sbps_all.size();++i)
{
double x=(radii_all[i]+radii_all[i-1])/2;
double ym=mv[i-1];
@@ -403,8 +403,8 @@ int main(int argc,char* argv[])
}
ofs_sbp<<"no no no"<<endl;
//bkg level
- double bkg_level=abs(f.get_param_value("bkg"));
- for(int i=0;i<sbps_all.size();++i)
+ bkg_level=abs(f.get_param_value("bkg"));
+ for(size_t i=0;i<sbps_all.size();++i)
{
double x=(radii_all[i]+radii_all[i-1])/2;
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
@@ -420,23 +420,23 @@ int main(int argc,char* argv[])
}
//resid
ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
+ //double y=sbps[i-1];
+ //double ye=sbpe[i-1];
double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
}
//zero level of resid
ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
- double ym=mv[i-1];
+ //double y=sbps[i-1];
+ //double ye=sbpe[i-1];
+ //double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl;
}
@@ -455,16 +455,16 @@ int main(int argc,char* argv[])
}
*/
- double lower,upper;
+ //double lower,upper;
double dr=1;
//calculate the mass profile
- const double G=6.673E-8;//cm^3 g^-1 s^-2
+ //const double G=6.673E-8;//cm^3 g^-1 s^-2
// Molecular weight per electron
// Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below
static const double mu=1.155;
static const double mp=1.67262158E-24;//g
static const double M_sun=1.98892E33;//g
- static const double k=1.38E-16;
+ //static const double k=1.38E-16;
ofstream ofs_mass("mass_int.qdp");
ofstream ofs_mass_dat("mass_int.dat");
@@ -494,14 +494,14 @@ int main(int argc,char* argv[])
double T_keV=Tprof(r);
double T1_keV=Tprof(r1);
- double T_K=T_keV*11604505.9;
- double T1_K=T1_keV*11604505.9;
+ //double T_K=T_keV*11604505.9;
+ //double T1_K=T1_keV*11604505.9;
double dlnT=log(T1_keV/T_keV);
double dlnr=log(r+dr)-log(r);
double dlnn=log(ne1/ne);
- double r_kpc=r_cm/kpc;
+ //double r_kpc=r_cm/kpc;
double r_Mpc=r_cm/Mpc;
//double M=-r_cm*T_K*k/G/mu/mp*(dlnT/dlnr+dlnn/dlnr);
//ref:http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W
diff --git a/mass_profile/fit_dbeta_sbp.cpp b/mass_profile/fit_dbeta_sbp.cpp
index 833059c..7fb1c7d 100644
--- a/mass_profile/fit_dbeta_sbp.cpp
+++ b/mass_profile/fit_dbeta_sbp.cpp
@@ -180,7 +180,7 @@ int main(int argc,char* argv[])
for(ifstream ifs(cfg.cfunc_file.c_str());;)
{
assert(ifs.is_open());
- double x,y,y1,y2;
+ double x,y;
ifs>>x>>y;
if(!ifs.good())
{
@@ -270,7 +270,7 @@ int main(int argc,char* argv[])
double n02=0;
double rc2=0;
double beta=0;
- double bkg=0;
+ double bkg_level=0;
if(tie_beta)
{
f.set_param_value("beta",.7);
@@ -362,7 +362,7 @@ int main(int argc,char* argv[])
//output the params
ofstream param_output("dbeta_param.txt");
//output the datasets and fitting results
- for(int i=0;i<f.get_num_params();++i)
+ for(size_t i=0;i<f.get_num_params();++i)
{
if(f.get_param_info(i).get_name()=="rc1")
{
@@ -423,27 +423,27 @@ int main(int argc,char* argv[])
ofs_sbp<<"log x"<<endl;
ofs_sbp<<"log y off"<<endl;
ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
double y=sbps[i-1];
double ye=sbpe[i-1];
- double ym=mv[i-1];
+ //double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<endl;
}
ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
+ //double y=sbps[i-1];
+ //double ye=sbpe[i-1];
double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<0<<endl;
}
//bkg
ofs_sbp<<"no no no"<<endl;
- double bkg_level=abs(f.get_param_value("bkg"));
- for(int i=0;i<sbps.size();++i)
+ bkg_level=abs(f.get_param_value("bkg"));
+ for(size_t i=0;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
@@ -466,22 +466,22 @@ int main(int argc,char* argv[])
}
//resid
ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
+ //double y=sbps[i-1];
+ //double ye=sbpe[i-1];
double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
}
//zero level in resid map
ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
- double ym=mv[i-1];
+ //double y=sbps[i-1];
+ //double ye=sbpe[i-1];
+ //double ym=mv[i-1];
ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl;
}
@@ -501,16 +501,16 @@ int main(int argc,char* argv[])
}
*/
- double lower,upper;
+ //double lower,upper;
double dr=1;
//calculate the mass profile
- const double G=6.673E-8;//cm^3 g^-1 s^-2
+ //const double G=6.673E-8;//cm^3 g^-1 s^-2
// Molecular weight per electron
// Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below
static const double mu=1.155;
static const double mp=1.67262158E-24;//g
static const double M_sun=1.98892E33;//g
- static const double k=1.38E-16;
+ //static const double k=1.38E-16;
ofstream ofs_mass("mass_int.qdp");
ofstream ofs_mass_dat("mass_int.dat");
@@ -536,26 +536,23 @@ int main(int argc,char* argv[])
gas_mass+=dmgas;
ofs_gas_mass<<r*cm_per_pixel/kpc<<"\t"<<gas_mass<<endl;
- double ne_beta1=dbeta_func(r,n01,rc1,beta1,
- 0,rc2,beta2);
+ double ne_beta1=dbeta_func(r,n01,rc1,beta1, 0,rc2,beta2);
- double ne_beta2=dbeta_func(r,0,rc1,beta1,
- n02,rc2,beta2);
+ double ne_beta2=dbeta_func(r,0,rc1,beta1, n02,rc2,beta2);
- double ne1=dbeta_func(r1,n01,rc1,beta1,
- n02,rc2,beta2);//cm^3
+ double ne1=dbeta_func(r1,n01,rc1,beta1, n02,rc2,beta2);//cm^3
double T_keV=Tprof(r);
double T1_keV=Tprof(r1);
- double T_K=T_keV*11604505.9;
- double T1_K=T1_keV*11604505.9;
+ //double T_K=T_keV*11604505.9;
+ //double T1_K=T1_keV*11604505.9;
double dlnT=log(T1_keV/T_keV);
double dlnr=log(r+dr)-log(r);
double dlnn=log(ne1/ne);
- double r_kpc=r_cm/kpc;
+ //double r_kpc=r_cm/kpc;
double r_Mpc=r_cm/Mpc;
//double M=-r_cm*T_K*k/G/mu/mp*(dlnT/dlnr+dlnn/dlnr);
//ref:http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W
diff --git a/mass_profile/fit_direct_beta.cpp b/mass_profile/fit_direct_beta.cpp
index cc1d7c9..ce9c63c 100644
--- a/mass_profile/fit_direct_beta.cpp
+++ b/mass_profile/fit_direct_beta.cpp
@@ -22,7 +22,7 @@ int main(int argc,char* argv[])
}
fitter<double,double,vector<double>,double,string> f;
-
+
f.set_statistic(chisq<double,double,vector<double>,double,string>());
f.set_opt_method(powell_method<double,vector<double> >());
f.set_model(beta1d<double>());
@@ -31,11 +31,11 @@ int main(int argc,char* argv[])
ifs>>dl;
f.load_data(dl.get_data_set());
f.fit();
-
+
double rmin=f.get_data_set().get_data(0).get_x();
double rmax=f.get_data_set().get_data(f.get_data_set().size()-1).get_x();
cout<<"read terr 1 2\nskip single\n";
- for(int i=0;i<f.get_data_set().size();++i)
+ for(size_t i=0;i<f.get_data_set().size();++i)
{
cout<<f.get_data_set().get_data(i).get_x()<<"\t"<<
-abs(f.get_data_set().get_data(i).get_x_lower_err())<<"\t"<<
@@ -44,7 +44,7 @@ int main(int argc,char* argv[])
-abs(f.get_data_set().get_data(i).get_y_lower_err())<<"\t"<<
abs(f.get_data_set().get_data(i).get_y_upper_err())<<endl;
-
+
}
cout<<"no no no\n";
@@ -53,7 +53,7 @@ int main(int argc,char* argv[])
cout<<i<<"\t0\t0\t"<<f.eval_model(i,f.get_all_params())<<"\t0\t0"<<endl;
}
- for(int i=0;i<f.get_num_params();++i)
+ for(size_t i=0;i<f.get_num_params();++i)
{
cerr<<f.get_param_info(i).get_name()<<"="<<
f.get_param_info(i).get_value()<<endl;
diff --git a/mass_profile/fit_lt_bpl.cpp b/mass_profile/fit_lt_bpl.cpp
index 5637dca..0c06b22 100644
--- a/mass_profile/fit_lt_bpl.cpp
+++ b/mass_profile/fit_lt_bpl.cpp
@@ -45,7 +45,7 @@ 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));
@@ -108,7 +108,7 @@ int main(int argc,char* argv[])
}
line+=" ";
istringstream iss(line);
-
+
if(line[0]=='#')
{
if(!is_first_nonono)
@@ -190,16 +190,16 @@ int main(int argc,char* argv[])
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>());
@@ -208,7 +208,7 @@ int main(int argc,char* argv[])
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);
@@ -225,7 +225,7 @@ int main(int argc,char* argv[])
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;
@@ -238,7 +238,7 @@ int main(int argc,char* argv[])
double sxl=0;
double syl=0;
double sxyl=0;
-
+
double sxxu=0;
double s1u=0;
double sxu=0;
@@ -246,7 +246,7 @@ int main(int argc,char* argv[])
double sxyu=0;
opt_utilities::default_data_set<double,double> ds1;
- for(int i=0;i<ds.size();++i)
+ 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();
@@ -254,14 +254,14 @@ int main(int argc,char* argv[])
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,
@@ -292,40 +292,40 @@ int main(int argc,char* argv[])
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(int i=0;i<p.size();++i)
+ 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(int i=0;i<mean_p.size();++i)
+ for(size_t i=0;i<mean_p.size();++i)
{
mean_p[i]/=cnt;
mean_p2[i]/=cnt;
diff --git a/mass_profile/fit_lt_pl.cpp b/mass_profile/fit_lt_pl.cpp
index e1686e5..fa37249 100644
--- a/mass_profile/fit_lt_pl.cpp
+++ b/mass_profile/fit_lt_pl.cpp
@@ -28,13 +28,13 @@ 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;
@@ -65,7 +65,7 @@ int main(int argc,char* argv[])
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]));
+ 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");
@@ -91,15 +91,15 @@ int main(int argc,char* argv[])
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)
@@ -163,7 +163,7 @@ int main(int argc,char* argv[])
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> >());
@@ -185,7 +185,7 @@ int main(int argc,char* argv[])
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;
@@ -196,7 +196,7 @@ int main(int argc,char* argv[])
++cnt;
cerr<<".";
opt_utilities::default_data_set<double,double> ds1;
- for(int i=0;i<ds.size();++i)
+ 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(),
@@ -212,7 +212,7 @@ int main(int argc,char* argv[])
//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");
@@ -234,6 +234,6 @@ int main(int argc,char* argv[])
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;
+ std::cout<<"Num of sources:"<<ds.size()<<endl;
}
diff --git a/mass_profile/fit_mt_bpl.cpp b/mass_profile/fit_mt_bpl.cpp
index 7cf22b2..e956f4f 100644
--- a/mass_profile/fit_mt_bpl.cpp
+++ b/mass_profile/fit_mt_bpl.cpp
@@ -45,7 +45,7 @@ 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));
@@ -107,7 +107,7 @@ int main(int argc,char* argv[])
}
line+=" ";
istringstream iss(line);
-
+
if(line[0]=='#')
{
if(!is_first_nonono)
@@ -145,7 +145,7 @@ int main(int argc,char* argv[])
cerr<<line<<endl;
continue;
}
-#endif
+#endif
if(std::abs(Mu)+std::abs(Ml)<M*.1)
{
@@ -204,16 +204,16 @@ int main(int argc,char* argv[])
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>());
@@ -222,7 +222,7 @@ int main(int argc,char* argv[])
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);
@@ -239,7 +239,7 @@ int main(int argc,char* argv[])
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;
@@ -252,7 +252,7 @@ int main(int argc,char* argv[])
double sxl=0;
double syl=0;
double sxyl=0;
-
+
double sxxu=0;
double s1u=0;
double sxu=0;
@@ -260,7 +260,7 @@ int main(int argc,char* argv[])
double sxyu=0;
opt_utilities::default_data_set<double,double> ds1;
- for(int i=0;i<ds.size();++i)
+ 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();
@@ -268,14 +268,14 @@ int main(int argc,char* argv[])
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,
@@ -306,40 +306,40 @@ int main(int argc,char* argv[])
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(int i=0;i<p.size();++i)
+ 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(int i=0;i<mean_p.size();++i)
+ for(size_t i=0;i<mean_p.size();++i)
{
mean_p[i]/=cnt;
mean_p2[i]/=cnt;
diff --git a/mass_profile/fit_mt_pl.cpp b/mass_profile/fit_mt_pl.cpp
index a2246a7..bfccaf4 100644
--- a/mass_profile/fit_mt_pl.cpp
+++ b/mass_profile/fit_mt_pl.cpp
@@ -92,7 +92,7 @@ int main(int argc,char* argv[])
}
line+=" ";
istringstream iss(line);
-
+
if(line[0]=='#')
{
if(!is_first_nonono)
@@ -129,7 +129,7 @@ int main(int argc,char* argv[])
cerr<<line<<endl;
continue;
}
-#endif
+#endif
if(std::abs(Mu)+std::abs(Ml)<M*.1)
{
double k=M*.1/(std::abs(Mu)+std::abs(Ml));
@@ -170,7 +170,7 @@ int main(int argc,char* argv[])
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> >());
@@ -203,12 +203,12 @@ int main(int argc,char* argv[])
ofs_resid<<"skip single"<<endl;
ofs_resid<<"ma 3 on 1"<<endl;
ofs_resid<<"log x"<<endl;
- for(int i=0;i<ds.size();++i)
+ 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 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;
@@ -223,7 +223,7 @@ int main(int argc,char* argv[])
++cnt;
cerr<<".";
opt_utilities::default_data_set<double,double> ds1;
- for(int i=0;i<ds.size();++i)
+ 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(),
@@ -238,7 +238,7 @@ int main(int argc,char* argv[])
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");
@@ -260,6 +260,6 @@ int main(int argc,char* argv[])
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;
+ std::cout<<"Num of sources:"<<ds.size()<<endl;
}
diff --git a/mass_profile/fit_nfw_mass.cpp b/mass_profile/fit_nfw_mass.cpp
index 48064be..fae74ef 100644
--- a/mass_profile/fit_nfw_mass.cpp
+++ b/mass_profile/fit_nfw_mass.cpp
@@ -1,7 +1,7 @@
/*
Fitting nfw mass profile model
Author: Junhua Gu
- Last modification 20120721
+ Last modification 20120721
*/
#include "nfw.hpp"
@@ -26,7 +26,7 @@ static double calc_critical_density(double z,
{
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;
+ const double H=H0*E;
return 3*H*H/8/pi/G;
}
@@ -53,14 +53,14 @@ int main(int argc,char* argv[])
ifstream ifs(argv[1]);
//cout<<"read serr 2"<<endl;
ofstream ofs_fit_result("nfw_fit_result.qdp");
-
+
ofs_fit_result<<"read serr 1 2"<<endl;
ofs_fit_result<<"skip single"<<endl;
ofs_fit_result<<"line off"<<endl;
ofs_fit_result<<"li on 2"<<endl;
ofs_fit_result<<"li on 4"<<endl;
ofs_fit_result<<"ls 2 on 4"<<endl;
-
+
ofs_fit_result<<"win 1"<<endl;
ofs_fit_result<<"yplot 1 2"<<endl;
ofs_fit_result<<"loc 0 0 1 1"<<endl;
@@ -108,7 +108,7 @@ int main(int argc,char* argv[])
vector<double> p=fit.fit();
//output parameters
ofstream ofs_param("nfw_param.txt");
- for(int i=0;i<fit.get_num_params();++i)
+ for(size_t i=0;i<fit.get_num_params();++i)
{
cout<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<endl;
ofs_param<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<endl;
@@ -117,11 +117,11 @@ int main(int argc,char* argv[])
ofs_param<<"reduced chi^2="<<fit.get_statistic_value()<<endl;
ofstream ofs_model("nfw_dump.qdp");
ofstream ofs_overdensity("overdensity.qdp");
- const double G=6.673E-8;//cm^3 g^-1 s^-2
- static const double mu=1.4074;
- static const double mp=1.67262158E-24;//g
+ //const double G=6.673E-8;//cm^3 g^-1 s^-2
+ //static const double mu=1.4074;
+ //static const double mp=1.67262158E-24;//g
static const double M_sun=1.98892E33;//g
- static const double k=1.38E-16;
+ //static const double k=1.38E-16;
double xmax=0;
for(double x=std::max(rmin_kpc,ds.get_data(0).get_x());;x+=1)
@@ -141,7 +141,7 @@ int main(int argc,char* argv[])
}
}
ofs_fit_result<<"no no no"<<endl;
- for(int i=0;i<ds.size();++i)
+ for(size_t i=0;i<ds.size();++i)
{
data<double,double> d=ds.get_data(i);
double x=d.get_x();
@@ -155,5 +155,5 @@ int main(int argc,char* argv[])
{
ofs_fit_result<<x<<"\t"<<0<<"\t"<<0<<"\t"<<0<<endl;
}
-
+
}
diff --git a/mass_profile/fit_nfw_sbp.cpp b/mass_profile/fit_nfw_sbp.cpp
index 6343dea..a7d0537 100644
--- a/mass_profile/fit_nfw_sbp.cpp
+++ b/mass_profile/fit_nfw_sbp.cpp
@@ -46,7 +46,7 @@ public:
*/
return spl.get_value(x);
}
-
+
//we need this function, when this object is performing a clone of itself
spline_func_obj* do_clone()const
{
@@ -59,7 +59,7 @@ public:
{
spl.push_point(x,y);
}
-
+
//before getting the intepolated value, the spline should be initialzied by calling this function
void gen_spline()
{
@@ -137,9 +137,9 @@ int main(int argc,char* argv[])
this will be convenient to calculate the volume of each spherical shell,
and can naturally ensure the annuli are adjacent with each other, with out any gaps.
-
+
*/
-
+
std::vector<double> radii;//to store radius
std::vector<double> sbps;//to store the surface brightness value
std::vector<double> sbpe;//to store the sbp error
@@ -148,8 +148,8 @@ int main(int argc,char* argv[])
About the format of the radius file:
the radius file contains only radius, separated by space or line feed (i.e., the <ENTER> key).
the unit should be pixel
-
- The number of radius can be larger than the number of annuli+1, the exceeded radius can be used
+
+ The number of radius can be larger than the number of annuli+1, the exceeded radius can be used
to calculate the influence of outer shells.
*/
int ncut=0;
@@ -208,7 +208,7 @@ int main(int argc,char* argv[])
for(ifstream ifs(arg_map["cfunc_file"].c_str());;)
{
assert(ifs.is_open());
- double x,y,y1,y2;
+ double x,y;
ifs>>x>>y;
if(!ifs.good())
{
@@ -250,7 +250,7 @@ int main(int argc,char* argv[])
nfw.set_cm_per_pixel(cm_per_pixel);
//define the temperature profile model
spline_func_obj tf;
-
+
for(ifstream ifs_tfunc(arg_map["T_file"].c_str());;)
{
assert(ifs_tfunc.is_open());
@@ -286,8 +286,8 @@ int main(int argc,char* argv[])
f.set_param_value("n0",n0);
f.set_param_value("rho0",rho0);
f.set_param_value("rs",rs);
-
-
+
+
f.set_param_value("bkg",bkg);
cout<<f.get_data_set().size()<<endl;
@@ -302,8 +302,8 @@ int main(int argc,char* argv[])
#endif
//output the parameters
ofstream param_output("nfw_param.txt");
-
- for(int i=0;i<f.get_num_params();++i)
+
+ for(size_t i=0;i<f.get_num_params();++i)
{
cout<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
@@ -326,7 +326,7 @@ int main(int argc,char* argv[])
//output the surface brightness profile
ofs_sbp<<"read serr 2"<<endl;
ofs_sbp<<"skip single"<<endl;
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
double y=sbps[i-1];
@@ -338,7 +338,7 @@ int main(int argc,char* argv[])
//output the electron density
mv=nfw.eval(radii,p);
ofstream ofs_rho("rho_fit.qdp");
- for(int i=1;i<sbps.size();++i)
+ for(size_t i=1;i<sbps.size();++i)
{
double x=(radii[i]+radii[i-1])/2;
double ym=mv[i-1];
@@ -354,7 +354,7 @@ int main(int argc,char* argv[])
}
//calculate the overdensity profile
ofstream ofs_overdensity("overdensity.qdp");
-
+
std::vector<double> radius_list;
std::vector<double> delta_list;
@@ -380,8 +380,8 @@ int main(int argc,char* argv[])
*/
ofs_overdensity<<r/kpc<<"\t"<<delta<<endl;
}
-
- for(int i=0;i<radius_list.size()-1;++i)
+
+ for(size_t i=0;i<radius_list.size()-1;++i)
{
double r=radius_list[i];
if(delta_list[i]>=200&&delta_list[i+1]<200)
diff --git a/mass_profile/fit_wang2012_model.cpp b/mass_profile/fit_wang2012_model.cpp
index 3fe14b6..f15a2b1 100644
--- a/mass_profile/fit_wang2012_model.cpp
+++ b/mass_profile/fit_wang2012_model.cpp
@@ -2,7 +2,7 @@
Fitting Jy Wang's temperature profile model
Author: Jingying Wang
Last modification 20120819
-
+
*/
#include "wang2012_model.hpp"
@@ -34,7 +34,7 @@ int main(int argc,char* argv[])
{
cm_per_pixel=atof(argv[3]);
}
-
+
//define the fitter
fitter<double,double,vector<double>,double,std::string> fit;
//define the data set
@@ -88,7 +88,7 @@ int main(int argc,char* argv[])
fit.set_statistic(chisq_object);
//fit.set_statistic(chisq<double,double,vector<double>,double,std::string>());
fit.set_model(wang2012_model<double>());
-
+
if(argc>=3&&std::string(argv[2])!="NONE")
{
std::vector<std::string> freeze_list;
@@ -123,7 +123,7 @@ int main(int argc,char* argv[])
{
freeze_param<double,double,std::vector<double>,std::string> fp(freeze_list[0]);
fit.set_param_modifier(fp);
- for(int i=1;i<freeze_list.size();++i)
+ for(size_t i=1;i<freeze_list.size();++i)
{
dynamic_cast<freeze_param<double,double,std::vector<double>,std::string>&>(fit.get_param_modifier())+=freeze_param<double,double,std::vector<double>,std::string>(freeze_list[i]);
}
@@ -147,11 +147,11 @@ int main(int argc,char* argv[])
}
#endif
//output parameters
- for(int i=0;i<fit.get_num_params();++i)
+ for(size_t i=0;i<fit.get_num_params();++i)
{
std::string pname=fit.get_param_info(i).get_name();
std::string pstatus=fit.get_model().report_param_status(pname);
-
+
if(pstatus==""||pstatus=="thawed")
{
pstatus="T";
@@ -168,20 +168,20 @@ int main(int argc,char* argv[])
}
#endif
}
-
+
//cout<<"T0"<<"\t"<<fit.get_param_value("T0")<<endl;
//cout<<"T1"<<"\t"<<fit.get_param_value("T1")<<endl;
//cout<<"xt"<<"\t"<<fit.get_param_value("xt")<<endl;
//cout<<"eta"<<"\t"<<fit.get_param_value("eta")<<endl;
//dump the data for checking
ofstream ofs_model("wang2012_dump.qdp");
-
-
+
+
min_r=0;
for(double x=min_r;x<3000;x+=10)
{
double model_value=fit.eval_model_raw(x,p);
-
+
ofs_model<<x<<"\t"<<model_value<<endl;
if(cm_per_pixel>0)
{
diff --git a/mass_profile/nfw_ne.hpp b/mass_profile/nfw_ne.hpp
index 5842bb8..8f23820 100644
--- a/mass_profile/nfw_ne.hpp
+++ b/mass_profile/nfw_ne.hpp
@@ -36,7 +36,7 @@ namespace opt_utilities
{
return rho0;
}
-
+
return nfw_mass_enclosed(r,rho0,rs)/(4.*pi/3*r*r*r);
}
@@ -46,7 +46,7 @@ namespace opt_utilities
const double Omega_m=.27)
{
const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m);
- const double H=H0*E;
+ const double H=H0*E;
return 3*H*H/8/pi/G;
}
@@ -66,7 +66,7 @@ namespace opt_utilities
nfw_ne()
:pTfunc(0),cm_per_pixel(1)
{
-
+
this->push_param_info(param_info<std::vector<T>,std::string>("rho0",1));//in mp
this->push_param_info(param_info<std::vector<T>,std::string>("rs",100));
this->push_param_info(param_info<std::vector<T>,std::string>("n0",.01));
@@ -89,7 +89,7 @@ namespace opt_utilities
this->push_param_info(param_info<std::vector<T>,std::string>("rs",rhs.get_param_info("rs").get_value()));
this->push_param_info(param_info<std::vector<T>,std::string>("n0",rhs.get_param_info("n0").get_value()));
}
-
+
//assignment operator
nfw_ne& operator=(const nfw_ne& rhs)
{
@@ -146,19 +146,19 @@ namespace opt_utilities
*/
std::vector<T> do_eval(const std::vector<T> & r,
const std::vector<T>& p)
- {
+ {
assert(pTfunc);
//const T kT_erg=k*5;
T rho0=std::abs(p[0])*mp;
T rs=std::abs(p[1]);
T n0=std::abs(p[2]);
T rs_cm=rs*cm_per_pixel;
-
+
std::vector<T> yvec(r.size());
const T kT_erg0=pTfunc->eval((r.at(0)+r.at(1))/2)*k;
//calculate the integration
#pragma omp parallel for schedule(dynamic)
- for(int i=0;i<r.size();++i)
+ for(size_t i=0;i<r.size();++i)
{
T r_cm=r[i]*cm_per_pixel;
T kT_erg=pTfunc->eval(r[i])*k;
@@ -170,10 +170,10 @@ namespace opt_utilities
//std::cout<<r_cm/1e20<<"\t"<<nfw_mass_enclosed(r_cm,rho0,rs_cm)/1e45<<std::endl;
//std::cout<<r_cm/1e20<<"\t"<<G*nfw_mass_enclosed(r_cm,rho0,rs_cm)*mu*mp/kT_erg/r_cm/r_cm<<std::endl;
}
-
+
std::vector<T> ydxvec(r.size()-1);
#pragma omp parallel for schedule(dynamic)
- for(int i=1;i<r.size();++i)
+ for(size_t i=1;i<r.size();++i)
{
T dr=r[i]-r[i-1];
T dr_cm=dr*cm_per_pixel;
@@ -183,7 +183,7 @@ namespace opt_utilities
//construct the result
std::vector<T> result(r.size()-1);
#pragma omp parallel for schedule(dynamic)
- for(int i=0;i<r.size()-1;++i)
+ for(size_t i=0;i<r.size()-1;++i)
{
T y=-ydxvec.at(i);
T kT_erg=pTfunc->eval(r[i])*k;
diff --git a/mass_profile/projector.hpp b/mass_profile/projector.hpp
index a9d5f2a..6f2f192 100644
--- a/mass_profile/projector.hpp
+++ b/mass_profile/projector.hpp
@@ -99,7 +99,7 @@ namespace opt_utilities
void attach_model(const model<std::vector<T>,std::vector<T>,std::vector<T> >& m)
{
this->clear_param_info();
- for(int i=0;i<m.get_num_params();++i)
+ for(size_t i=0;i<m.get_num_params();++i)
{
this->push_param_info(m.get_param_info(i));
}
@@ -167,7 +167,7 @@ namespace opt_utilities
}
}
std::vector<T> p2(p1.size()-1);
- for(int i=0;i<p1.size()-1;++i)
+ for(size_t i=0;i<p1.size()-1;++i)
{
p2.at(i)=p1[i];
}
@@ -183,10 +183,10 @@ namespace opt_utilities
std::vector<T> unprojected(pmodel->eval(x,p));
std::vector<T> projected(unprojected.size());
- for(int nrad=0;nrad<x.size()-1;++nrad)
+ for(size_t nrad=0;nrad<x.size()-1;++nrad)
{
double v1=0;
- for(int nsph=nrad;nsph<x.size()-1;++nsph)
+ for(size_t nsph=nrad;nsph<x.size()-1;++nsph)
{
double v=calc_v(x,nsph,nrad)*cm_per_pixel*cm_per_pixel*cm_per_pixel;
v1+=v;
diff --git a/mass_profile/vchisq.hpp b/mass_profile/vchisq.hpp
index 0780ce8..c3d38b8 100644
--- a/mass_profile/vchisq.hpp
+++ b/mass_profile/vchisq.hpp
@@ -28,7 +28,7 @@ namespace opt_utilities
int n;
bool limit_bound;
typedef std::vector<T> Tp;
-
+
vchisq<T>* do_clone()const
{
return new vchisq<T>(*this);
@@ -38,7 +38,7 @@ namespace opt_utilities
{
return "chi^2 statistic";
}
-
+
public:
void verbose(bool v)
{
@@ -49,7 +49,7 @@ namespace opt_utilities
{
limit_bound=true;
}
-
+
void clear_limit()
{
limit_bound=false;
@@ -58,8 +58,8 @@ namespace opt_utilities
vchisq()
:verb(false),limit_bound(false)
{}
-
-
+
+
T do_eval(const std::vector<T>& p)
{
@@ -71,7 +71,7 @@ namespace opt_utilities
}
}
T result(0);
-
+
std::vector<float> vx;
std::vector<float> vy;
std::vector<float> vye;
@@ -80,7 +80,7 @@ namespace opt_utilities
if(verb)
{
n++;
-
+
if(n%100==0)
{
vx.resize(this->get_data_set().get_data(0).get_y().size());
@@ -88,24 +88,24 @@ namespace opt_utilities
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<double> y_model(this->eval_model(this->get_data_set().get_data(i).get_x(),p));
const std::vector<double>& y=this->get_data_set().get_data(i).get_y();
const std::vector<double>& ye=this->get_data_set().get_data(i).get_y_lower_err();
- for(int j=0;j<y.size();++j)
+ for(size_t j=0;j<y.size();++j)
{
double chi=(y_model[j]-y[j])/ye[j];
result+=chi*chi;
}
-
-
+
+
if(verb&&n%100==0)
{
-
- for(int j=0;j<y.size();++j)
+
+ for(size_t j=0;j<y.size();++j)
{
vx.at(j)=((this->get_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]);
@@ -121,9 +121,9 @@ namespace opt_utilities
my[j]=log10(my[j]);
}
-
+
}
-
+
}
if(verb)
{
@@ -148,7 +148,7 @@ namespace opt_utilities
pr.plot_line(vx,my);
}
}
-
+
return result;
}
};
Copyright © 2017-2018 Aaron LI. All Rights Reserved.