noe/pipeline
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Remove log in Brho, clean and comment

Browse Source
This commit is contained in:
Noe Brucy
2020-06-23 10:31:18 +02:00
committed by Noe Brucy
parent d92fb888f9
commit 850cc4991e
2 changed files with 94 additions and 54 deletions
Download Patch File Download Diff File Expand all files Collapse all files
plotter.py
+55 -3
View File
@@ -142,7 +142,9 @@ class Plotter(Aggregator, BaseProcessor):
self.save = None self.save = None
def _not_self_dep(self, name, dep, dep_arg, overwrite, **kwargs): def _not_self_dep(self, name, dep, dep_arg, overwrite, **kwargs):
"""""" """
Check if the dependency belongs to the plotter object or to another one (comp, pp, ..)
"""
if dep in self.comp.rules: if dep in self.comp.rules:
done = self.comp.process(dep, dep_arg, overwrite, overwrite) done = self.comp.process(dep, dep_arg, overwrite, overwrite)
self.just_done.extend(done) self.just_done.extend(done)
@@ -150,6 +152,9 @@ class Plotter(Aggregator, BaseProcessor):
super(Plotter, self)._not_self_dep(name, dep, dep_arg, overwrite, **kwargs) super(Plotter, self)._not_self_dep(name, dep, dep_arg, overwrite, **kwargs)
def _needs_computation(self, overwrite, plot_filename): def _needs_computation(self, overwrite, plot_filename):
"""
Returns true if the plot needs to be redone
"""
return ( return (
self.pp_params.out.interactive self.pp_params.out.interactive
or overwrite or overwrite
@@ -167,6 +172,9 @@ class Plotter(Aggregator, BaseProcessor):
from_cells=False, from_cells=False,
**kwargs **kwargs
): ):
"""
Generic method to process a rule, with its dependencies
"""
if not arg is None: if not arg is None:
name_full = name + "_" + str(arg) name_full = name + "_" + str(arg)
else: else:
@@ -262,6 +270,9 @@ class Plotter(Aggregator, BaseProcessor):
self._log("{} is not valid in this context".format(name_full), "ERROR") self._log("{} is not valid in this context".format(name_full), "ERROR")
def _plot_rule(self, rule, save, arg, plot_filename, overwrite, ax, **kwargs): def _plot_rule(self, rule, save, arg, plot_filename, overwrite, ax, **kwargs):
"""
Once all dependencies are met, actually process the rule
"""
P.sca(ax) P.sca(ax)
if self._needs_computation(overwrite, plot_filename): if self._needs_computation(overwrite, plot_filename):
rule.plot(save, arg, **kwargs) rule.plot(save, arg, **kwargs)
@@ -278,7 +289,9 @@ class Plotter(Aggregator, BaseProcessor):
self._log("Plot {} is already done, skipping...".format(plot_filename)) self._log("Plot {} is already done, skipping...".format(plot_filename))
def _find_filename(self, name_full, run=None, num=None, fmt=None): def _find_filename(self, name_full, run=None, num=None, fmt=None):
"""
Determine a filename based on rule name, run, output and parameters
"""
tag_name = self.pp_params.out.tag tag_name = self.pp_params.out.tag
if fmt is None and self.pp_params.out.fmt == "": if fmt is None and self.pp_params.out.fmt == "":
@@ -309,6 +322,10 @@ class Plotter(Aggregator, BaseProcessor):
) )
def _label_run(self, run, node, label, nml_key): def _label_run(self, run, node, label, nml_key):
"""
Set up a label for the run from the namelist and parameters
"""
def get_label_nml(nml_key): def get_label_nml(nml_key):
prop_name = os.path.basename(nml_key) prop_name = os.path.basename(nml_key)
if prop_name in self.label_convert: if prop_name in self.label_convert:
@@ -344,6 +361,9 @@ class Plotter(Aggregator, BaseProcessor):
return label_run return label_run
def _ax_label_unit(self, node, label, unit, unit_coeff): def _ax_label_unit(self, node, label, unit, unit_coeff):
"""
Find appropriate labels for axis
"""
if label is None: if label is None:
if "label" in node._v_attrs: if "label" in node._v_attrs:
label = node._v_attrs.label label = node._v_attrs.label
@@ -393,6 +413,9 @@ class Plotter(Aggregator, BaseProcessor):
autoscale=True, autoscale=True,
**kwargs **kwargs
): ):
"""
Plot data on a map
"""
ax_h = self._axes_h[ax_los] ax_h = self._axes_h[ax_los]
ax_v = self._axes_v[ax_los] ax_v = self._axes_v[ax_los]
@@ -457,6 +480,9 @@ class Plotter(Aggregator, BaseProcessor):
plot_overlay(ax_los) plot_overlay(ax_los)
def _overlay_levels(self, ax_los): def _overlay_levels(self, ax_los):
"""
Add an overlay : AMR levels
"""
map_level = self.save.get_node("/maps/{}_{}".format("levels", ax_los)).read() map_level = self.save.get_node("/maps/{}_{}".format("levels", ax_los)).read()
# Computing linewidths # Computing linewidths
levels_ar = np.arange(np.min(map_level), np.max(map_level) + 1) levels_ar = np.arange(np.min(map_level), np.max(map_level) + 1)
@@ -482,6 +508,9 @@ class Plotter(Aggregator, BaseProcessor):
def _overlay_speed( def _overlay_speed(
self, ax_los, unit=cst.km_s, unit_coeff=1.0, key_v=None, **kwargs self, ax_los, unit=cst.km_s, unit_coeff=1.0, key_v=None, **kwargs
): ):
"""
Add an overlay : velocity vector field
"""
ax_h = self._axes_h[ax_los] ax_h = self._axes_h[ax_los]
ax_v = self._axes_v[ax_los] ax_v = self._axes_v[ax_los]
dmap_vh_node = self.save.get_node("/maps/speed_h_{}".format(ax_los)) dmap_vh_node = self.save.get_node("/maps/speed_h_{}".format(ax_los))
@@ -524,6 +553,9 @@ class Plotter(Aggregator, BaseProcessor):
) )
def _overlay_B(self, ax_los, **kwargs): def _overlay_B(self, ax_los, **kwargs):
"""
Add an overlay : magnetic streamlines
"""
ax_h = self._axes_h[ax_los] ax_h = self._axes_h[ax_los]
ax_v = self._axes_v[ax_los] ax_v = self._axes_v[ax_los]
dmap_Bh_node = self.save.get_node("/maps/B_h_{}".format(ax_los)) dmap_Bh_node = self.save.get_node("/maps/B_h_{}".format(ax_los))
@@ -552,6 +584,9 @@ class Plotter(Aggregator, BaseProcessor):
P.streamplot(hh, vv, map_Bh_red, map_Bv_red) P.streamplot(hh, vv, map_Bh_red, map_Bv_red)
def _plot_radial(self, name, ax_los, label=None, xlog=False, ylog=False): def _plot_radial(self, name, ax_los, label=None, xlog=False, ylog=False):
"""
Plot a radial profile (for disks, OUTDATED)
"""
radial_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read() radial_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read()
bin_centers = 0.5 * (radial_bins[1:] + radial_bins[:-1]) bin_centers = 0.5 * (radial_bins[1:] + radial_bins[:-1])
@@ -593,7 +628,9 @@ class Plotter(Aggregator, BaseProcessor):
fitlabel=None, fitlabel=None,
**kwargs **kwargs
): ):
"""
Plot an histogram (PDF, etc ...)
"""
if not ax_los is None: if not ax_los is None:
name = name + "_" + ax_los name = name + "_" + ax_los
@@ -707,6 +744,9 @@ class Plotter(Aggregator, BaseProcessor):
subname_y=None, subname_y=None,
**kwargs **kwargs
): ):
"""
Generic plot routine, with name_x and name_y two path in the hdf5 file
"""
if not node_arg is None: if not node_arg is None:
name_x, name_y = name_x + "_" + node_arg, name_y + "_" + node_arg name_x, name_y = name_x + "_" + node_arg, name_y + "_" + node_arg
@@ -856,12 +896,18 @@ class Plotter(Aggregator, BaseProcessor):
hdf5_y.close() hdf5_y.close()
def _pspec(self, name, **kwargs): def _pspec(self, name, **kwargs):
"""
Plot power spectrum
"""
del kwargs["run"] del kwargs["run"]
file_pspec = self.save.get_node("/hdf5/pspec").read() file_pspec = self.save.get_node("/hdf5/pspec").read()
num = self.save.root._v_attrs.num num = self.save.root._v_attrs.num
getattr(pspec_read, "pspec_" + name)(file_pspec, ".", num, **kwargs) getattr(pspec_read, "pspec_" + name)(file_pspec, ".", num, **kwargs)
def _overlay_fit(self, x, y, yerr=None, kind="linear", label=None, **kwargs): def _overlay_fit(self, x, y, yerr=None, kind="linear", label=None, **kwargs):
"""
Add an overlay : fit a curve, linear or powerlaw
"""
if kind == "linear": if kind == "linear":
if yerr is None: if yerr is None:
(a, b, rho, _map_rule, stderr) = linregress(x, y) (a, b, rho, _map_rule, stderr) = linregress(x, y)
@@ -917,6 +963,9 @@ class Plotter(Aggregator, BaseProcessor):
P.plot(x, (10 ** b) * x ** a, label=label, **kwargs) P.plot(x, (10 ** b) * x ** a, label=label, **kwargs)
def overlay_kennicutt(self, n0, step): def overlay_kennicutt(self, n0, step):
"""
Add an overlay : kennicutt mass accretion
"""
P.grid(False) P.grid(False)
ylim = P.ylim() ylim = P.ylim()
(tmin, tmax) = P.xlim() (tmin, tmax) = P.xlim()
@@ -934,6 +983,9 @@ class Plotter(Aggregator, BaseProcessor):
P.ylim(ylim) P.ylim(ylim)
def def_rules(self): def def_rules(self):
"""
That is where rules are defined
"""
self.rules = { self.rules = {
# Generic rules # Generic rules
"plot": PlotRule( "plot": PlotRule(
postprocessor.py
+39 -51
View File
@@ -8,7 +8,7 @@ vol_func = lambda dset: dset["dx"] ** 3 # Volume function
getter_T = lambda dset: dset["P"] / dset["rho"] # Temperature getter_T = lambda dset: dset["P"] / dset["rho"] # Temperature
def _getter_abs_cos_vB(dset): def getter_abs_cos_vB(dset):
B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1)) B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1))
v_norm = np.sqrt(np.sum(dset["vel"] ** 2, axis=1)) v_norm = np.sqrt(np.sum(dset["vel"] ** 2, axis=1))
# Compute the dot product in each cell # Compute the dot product in each cell
@@ -16,12 +16,12 @@ def _getter_abs_cos_vB(dset):
return np.abs(dot_prod) / (v_norm * B_norm) return np.abs(dot_prod) / (v_norm * B_norm)
def _getter_B_int(dset): def getter_B_int(dset):
B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1)) B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1))
return B_norm return B_norm
def _getter_rho(dset): def getter_rho(dset):
return dset["rho"] return dset["rho"]
@@ -302,27 +302,44 @@ class PostProcessor(HDF5Container):
centers = 0.5 * (edges[1:] + edges[:-1]) centers = 0.5 * (edges[1:] + edges[:-1])
return (np.stack([values, centers]), {"logbins": logbins}) return (np.stack([values, centers]), {"logbins": logbins})
def _Brho(self, bins=100): def _Brho(self, bins=100, logbins=True):
"""
Mean of B in rho bins
"""
self.load_cells() self.load_cells()
B = _getter_B_int(self.cells) B = getter_B_int(self.cells)
rho = _getter_rho(self.cells) rho = getter_rho(self.cells)
B = np.log10(B)
rho = np.log10(rho) if logbins:
abscisse = np.linspace(np.min(rho), np.max(rho), bins) rho_bins = np.logspace(
np.log10(np.min(rho)), np.log10(np.max(rho)), bins, base=10
)
else:
rho_bins = np.linspace(np.min(rho), np.max(rho), bins)
weights = mass_func(self.cells) weights = mass_func(self.cells)
# For each cell, bin_number contains the number of the bins it belongs to
bin_number = np.zeros(len(B)) bin_number = np.zeros(len(B))
for rho_min in abscisse[:-1]:
# Go through the min value of rho of each bin
for rho_min in rho_bins[:-1]:
bin_number = bin_number + (rho > rho_min).astype(int) bin_number = bin_number + (rho > rho_min).astype(int)
B_mean = np.zeros(len(abscisse) - 1) # Compute the mean in each bin
B_mean = np.zeros(len(rho_bins) - 1)
for i in range(len(B_mean)): for i in range(len(B_mean)):
B_mean[i] = np.mean(B[bin_number == i]) B_mean[i] = np.mean(B[bin_number == i])
values, edges = np.histogram(B, abscisse, weights=weights)
centers = 0.5 * (edges[1:] + edges[:-1]) # Get the center of each bin
print("centers") if logbins:
centers = 10 ** (0.5 * (np.log10(rho_bins[1:]) + np.log10(rho_bins[:-1])))
else:
centers = 0.5 * (rho_bins[1:] + rho_bins[:-1])
if self.pp_params.process.unload_cells: if self.pp_params.process.unload_cells:
self.unload_cells() self.unload_cells()
return {"rho": centers, "B": B_mean} return ({"rho": centers, "B": B_mean}, {"logbins": logbins})
def _mwa_sigma(self, axes=["x", "y", "z"]): def _mwa_sigma(self, axes=["x", "y", "z"]):
mw_speed = self.save.get_node("/globals/mwa_speed").read() mw_speed = self.save.get_node("/globals/mwa_speed").read()
@@ -381,6 +398,9 @@ class PostProcessor(HDF5Container):
return self._vector_v("Br", "unit_mag", ax_los, z) return self._vector_v("Br", "unit_mag", ax_los, z)
def _B_int(self, ax_los, z=0.0): def _B_int(self, ax_los, z=0.0):
"""
Slice ont the intensity of the magnétic field
"""
B_op = ScalarOperator( B_op = ScalarOperator(
lambda dset: np.sqrt(np.sum(dset["Br"] ** 2, axis=1)), lambda dset: np.sqrt(np.sum(dset["Br"] ** 2, axis=1)),
self._ro.info["unit_mag"], self._ro.info["unit_mag"],
@@ -606,38 +626,6 @@ class PostProcessor(HDF5Container):
pdf.attrs.var = np.var pdf.attrs.var = np.var
return True return True
# def kappa(self, mask):
# mask_flat = mask.flatten()
# fmap = self.get_value("/maps/fluct_coldens_z")
# nb_cells = np.sum(self.mask_flat)
# values, edges = np.histogram(np.log10(fmap.flatten()[mask_flat]),
# self.pp_params.pdf.nb_bin,
# weights = np.ones(nb_cells) / nb_cells)
# centers = 0.5 * (edges[1:] + edges[:-1])
# mask_fit = ((centers > self.pp.pp_params.pdf.xmin_fit) &
# (centers < self.pp.pp_params.pdf.xmax_fit) &
# (values > 0))
# (a, b, r, p, err) = linregress(centers[mask_fit], np.log10(values[mask_fit]))
# return a, err
# def gamma(self, mask):
# rho3d = np.log10(self.cells["rho"])
# P3d = np.log10(np.sqrt(self.cells["P"]))
# xy3d = self.pp.cells["pos"][:, :2] - 0.5
# self.rr3d = np.sqrt(np.sum((self.xy3d)**2, axis=1))
# nb_cells = np.sum(self.mask_flat)
# values, edges = np.histogram(np.log10(fmap.flatten()[mask_flat]),
# self.pp_params.pdf.nb_bin,
# weights = np.ones(nb_cells) / nb_cells)
# centers = 0.5 * (edges[1:] + edges[:-1])
# mask_fit = ((centers > self.pp.pp_params.pdf.xmin_fit) &
# (centers < self.pp.pp_params.pdf.xmax_fit) &
# (values > 0))
# (a, b, r, p, err) = linregress(centers[mask_fit], np.log10(values[mask_fit]))
# return a, err
def _sinks(self): def _sinks(self):
csv_name = ( csv_name = (
self.path self.path
@@ -829,17 +817,17 @@ class PostProcessor(HDF5Container):
), ),
"cos_pdf": Rule( "cos_pdf": Rule(
self, self,
partial(self._vol_pdf, _getter_abs_cos_vB, logbins=False), partial(self._vol_pdf, getter_abs_cos_vB, logbins=False),
"Global cos-PDF", "Global cos-PDF",
"/hist", "/hist",
unit=cst.none, unit=cst.none,
), ),
"Brho": Rule( "Brho": Rule(
self, self,
partial(self._Brho), self._Brho,
"Brho average", "Average of B as a function of rho",
"/datasets", "/datasets",
unit=self.info["unit_mag"], unit={"rho": self.info["unit_density"], "B": self.info["unit_mag"]},
), ),
# Profiles # Profiles
"axis": Rule( "axis": Rule(