summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorAaron LI <aaronly.me@outlook.com>2016-07-10 22:02:25 +0800
committerAaron LI <aaronly.me@outlook.com>2016-07-10 22:02:25 +0800
commit2f89d7da79e4d7984575d12d88d8d42b44d45d05 (patch)
tree9f988523a3589d4016bfbb3ba29955540099f843
parent5eee6efd5fb0dd5d9c72d12bbbaf73421d462339 (diff)
downloadcexcess-2f89d7da79e4d7984575d12d88d8d42b44d45d05.tar.bz2
deproject_sbp.py: use "SmoothSpline" from module "spline.py"
-rwxr-xr-xdeproject_sbp.py83
1 files changed, 18 insertions, 65 deletions
diff --git a/deproject_sbp.py b/deproject_sbp.py
index a11078d..dc8601d 100755
--- a/deproject_sbp.py
+++ b/deproject_sbp.py
@@ -2,9 +2,11 @@
#
# Aaron LI
# Created: 2016-06-10
-# Updated: 2016-07-04
+# Updated: 2016-07-10
#
# Change logs:
+# 2016-07-10:
+# * Use class 'SmoothSpline' from module
# 2016-07-04:
# * Use model's "report()" method
# * Add config "sbpexp_rcut"
@@ -176,12 +178,10 @@ from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from matplotlib.figure import Figure
from configobj import ConfigObj
-import rpy2.robjects as ro
-from rpy2.robjects.packages import importr
-
from astro_params import AstroParams, ChandraPixel
from projection import Projection
from fitting_models import PLCModel
+from spline import SmoothSpline
plt.style.use("ggplot")
@@ -436,12 +436,8 @@ class BrightnessProfile:
rho_gas = None
# fitted smoothing spline to the SBP
s_spline = None
- s_spline_log10 = None
# fitted smoothing spline to the cooling function profile
cf_spline = None
- cf_spline_log10 = None
- # call R through `rpy2`
- mgcv = importr("mgcv")
def __init__(self, sbp_data, cf_data, z):
self.load_data(data=sbp_data)
@@ -483,47 +479,22 @@ class BrightnessProfile:
def get_radius(self):
return (self.r.copy(), self.r_err.copy())
- def fit_spline(self, spline, log10=True):
- """
- Fit a smoothing line to the specified spline data,
- by calling 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 == "brightness":
- if log10:
- x = ro.FloatVector(np.log10(self.r))
- y = ro.FloatVector(np.log10(self.s))
- self.s_spline_log10 = True
- else:
- x = ro.FloatVector(self.r)
- y = ro.FloatVector(self.s)
- self.s_spline_log10 = False
- weights = ro.FloatVector(self.s / self.s_err)
- self.s_spline = self.mgcv.gam(
- ro.Formula("y ~ s(x, bs='ps')"),
- weights=weights,
- data=ro.DataFrame({"x": x, "y": y}),
- method="REML")
+ x = self.r
+ y = self.s
+ weights = self.s / self.s_err
+ spl = "s_spline"
elif spline == "cooling_function":
- if log10:
- x = ro.FloatVector(np.log10(self.cf_radius))
- y = ro.FloatVector(np.log10(self.cf_value))
- self.cf_spline_log10 = True
- else:
- x = ro.FloatVector(self.cf_radius)
- y = ro.FloatVector(self.cf_value)
- self.cf_spline_log10 = False
- self.cf_spline = self.mgcv.gam(
- ro.Formula("y ~ s(x, bs='ps')"),
- data=ro.DataFrame({"x": x, "y": y}),
- method="REML")
+ x = self.cf_radius
+ y = self.cf_value
+ weights = None
+ spl = "cf_spline"
+ setattr(self, spl, SmoothSpline(x=x, y=y, weights=weights))
+ getattr(self, spl).fit(log10=log10)
def eval_spline(self, spline, x):
"""
@@ -531,31 +502,13 @@ class BrightnessProfile:
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 == "brightness":
spl = self.s_spline
- log10 = self.s_spline_log10
elif spline == "cooling_function":
spl = self.cf_spline
- log10 = self.cf_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 calc_electron_density(self):
"""
@@ -566,9 +519,9 @@ class BrightnessProfile:
unit: [ cm^-3 ] if the units converted for input data
"""
if self.s_spline is None:
- self.fit_spline(spline="brightness", log10=True)
+ self.fit_spline(spline="brightness", log10=["x", "y"])
if self.cf_spline is None:
- self.fit_spline(spline="cooling_function", log10=False)
+ self.fit_spline(spline="cooling_function", log10=[])
#
s_new = self.eval_spline(spline="brightness", x=self.r)
cf_new = self.eval_spline(spline="cooling_function", x=self.r)