summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorAaron LI <aaronly.me@outlook.com>2016-07-13 17:11:06 +0800
committerAaron LI <aaronly.me@outlook.com>2016-07-13 17:11:06 +0800
commitf27f9326554802ef21250587f0f2252c6a1d7139 (patch)
treec48af54713b029049c410be30df1dc3b35036d28
parentbbe0a0cbfae52205601b3d1d858b26d2bfb9e32d (diff)
downloadcexcess-f27f9326554802ef21250587f0f2252c6a1d7139.tar.bz2
calc_overdensity.py: Use class "SmoothSpline" from module "spline.py"
-rwxr-xr-xcalc_overdensity.py78
1 files changed, 16 insertions, 62 deletions
diff --git a/calc_overdensity.py b/calc_overdensity.py
index 562311e..501f5dd 100755
--- a/calc_overdensity.py
+++ b/calc_overdensity.py
@@ -2,9 +2,11 @@
#
# Aaron LI
# Created: 2016-06-30
-# Updated: 2016-07-11
+# Updated: 2016-07-13
#
# Change logs:
+# 2016-07-13:
+# * Use class 'SmoothSpline' from module 'spline.py'
# 2016-07-11:
# * Use a default config to allow a minimal user config
# 2016-07-04:
@@ -35,10 +37,8 @@ import astropy.units as au
from astropy.cosmology import FlatLambdaCDM
from configobj import ConfigObj
-import rpy2.robjects as ro
-from rpy2.robjects.packages import importr
-
from astro_params import AstroParams
+from spline import SmoothSpline
config_default = """
@@ -83,11 +83,7 @@ class MassProfile:
redshift = None
# fitted smoothing spline
m_spline = None
- m_spline_log10 = None
od_spline = None
- od_spline_log10 = None
- # call R through `rpy2`
- mgcv = importr("mgcv")
def __init__(self, mass, mass_type="total"):
self.load_data(data=mass, mass_type=mass_type)
@@ -130,7 +126,7 @@ class MassProfile:
Calculate the radius at which the overdensity is delta.
"""
if self.od_spline is None:
- self.fit_spline(spline="overdensity", log10=True)
+ self.fit_spline(spline="overdensity", log10=["x", "y"])
if min(self.overdensity) > delta:
raise ValueError("min(overdensity) > %d" % delta)
r = optimize.newton(
@@ -140,7 +136,7 @@ class MassProfile:
def calc_mass_delta(self, r_delta):
if self.m_spline is None:
- self.fit_spline(spline="mass", log10=True)
+ self.fit_spline(spline="mass", log10=["x", "y"])
return self.eval_spline(spline="mass", x=r_delta)
def save(self, outfile):
@@ -153,49 +149,24 @@ class MassProfile:
header = "radius[kpc] radius_err[kpc] overdensity"
np.savetxt(outfile, data, header=header)
- def fit_spline(self, spline, log10):
- """
- Fit a smoothing line to the specified spline data,
- by utilizing the R `mgcv::gam()` function.
-
- If 'log10' is True, the input data are first applied the log-scale
- transform, and then fitted by the smoothing spline.
-
- The fitted spline allows extrapolation.
- """
+ def fit_spline(self, spline, log10=[]):
if spline not in self.SPLINES:
raise ValueError("invalid spline: %s" % spline)
#
if spline == "mass":
# input gas/total mass profile
- if log10:
- x = ro.FloatVector(np.log10(self.r))
- y = ro.FloatVector(np.log10(self.m))
- self.m_spline_log10 = True
- else:
- x = ro.FloatVector(self.r)
- y = ro.FloatVector(self.m)
- self.m_spline_log10 = False
- self.m_spline = self.mgcv.gam(
- ro.Formula("y ~ s(x, bs='ps')"),
- data=ro.DataFrame({"x": x, "y": y}),
- method="REML")
+ x = self.r
+ y = self.m
+ spl = "m_spline"
elif spline == "overdensity":
# calculated overdensity profile
- if log10:
- x = ro.FloatVector(np.log10(self.r))
- y = ro.FloatVector(np.log10(self.overdensity))
- self.od_spline_log10 = True
- else:
- x = ro.FloatVector(self.radius)
- y = ro.FloatVector(self.rho_total)
- self.od_spline_log10 = False
- self.od_spline = self.mgcv.gam(
- ro.Formula("y ~ s(x, bs='ps')"),
- data=ro.DataFrame({"x": x, "y": y}),
- method="REML")
+ x = self.radius
+ y = self.rho_total
+ spl = "od_spline"
else:
raise ValueError("invalid spline: %s" % spline)
+ setattr(self, spl, SmoothSpline(x=x, y=y))
+ getattr(self, spl).fit(log10=log10)
def eval_spline(self, spline, x):
"""
@@ -203,31 +174,14 @@ class MassProfile:
Also check whether the spline was fitted in the log-scale space,
and transform the evaluated values if necessary.
"""
- x = np.array(x)
- if x.shape == ():
- x = x.reshape((1,))
if spline == "mass":
spl = self.m_spline
- log10 = self.m_spline_log10
elif spline == "overdensity":
spl = self.od_spline
- log10 = self.od_spline_log10
else:
raise ValueError("invalid spline: %s" % spline)
#
- if log10:
- x_new = ro.ListVector({"x": ro.FloatVector(np.log10(x))})
- y_pred = self.mgcv.predict_gam(spl, newdata=x_new)
- y_pred = 10 ** np.array(y_pred)
- else:
- x_new = ro.ListVector({"x": ro.FloatVector(x)})
- y_pred = self.mgcv.predict_gam(spl, newdata=x_new)
- y_pred = np.array(y_pred)
- #
- if len(y_pred) == 1:
- return y_pred[0]
- else:
- return y_pred
+ return spl.eval(x)
def main():