diff --git a/plotter.py b/plotter.py index bc8e54c..2211e1c 100644 --- a/plotter.py +++ b/plotter.py @@ -142,7 +142,9 @@ class Plotter(Aggregator, BaseProcessor): self.save = None 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: done = self.comp.process(dep, dep_arg, overwrite, overwrite) 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) def _needs_computation(self, overwrite, plot_filename): + """ + Returns true if the plot needs to be redone + """ return ( self.pp_params.out.interactive or overwrite @@ -167,6 +172,9 @@ class Plotter(Aggregator, BaseProcessor): from_cells=False, **kwargs ): + """ + Generic method to process a rule, with its dependencies + """ if not arg is None: name_full = name + "_" + str(arg) else: @@ -262,6 +270,9 @@ class Plotter(Aggregator, BaseProcessor): self._log("{} is not valid in this context".format(name_full), "ERROR") def _plot_rule(self, rule, save, arg, plot_filename, overwrite, ax, **kwargs): + """ + Once all dependencies are met, actually process the rule + """ P.sca(ax) if self._needs_computation(overwrite, plot_filename): rule.plot(save, arg, **kwargs) @@ -278,7 +289,9 @@ class Plotter(Aggregator, BaseProcessor): self._log("Plot {} is already done, skipping...".format(plot_filename)) 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 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): + """ + Set up a label for the run from the namelist and parameters + """ + def get_label_nml(nml_key): prop_name = os.path.basename(nml_key) if prop_name in self.label_convert: @@ -344,6 +361,9 @@ class Plotter(Aggregator, BaseProcessor): return label_run def _ax_label_unit(self, node, label, unit, unit_coeff): + """ + Find appropriate labels for axis + """ if label is None: if "label" in node._v_attrs: label = node._v_attrs.label @@ -393,6 +413,9 @@ class Plotter(Aggregator, BaseProcessor): autoscale=True, **kwargs ): + """ + Plot data on a map + """ ax_h = self._axes_h[ax_los] ax_v = self._axes_v[ax_los] @@ -457,6 +480,9 @@ class Plotter(Aggregator, BaseProcessor): plot_overlay(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() # Computing linewidths levels_ar = np.arange(np.min(map_level), np.max(map_level) + 1) @@ -482,6 +508,9 @@ class Plotter(Aggregator, BaseProcessor): def _overlay_speed( 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_v = self._axes_v[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): + """ + Add an overlay : magnetic streamlines + """ ax_h = self._axes_h[ax_los] ax_v = self._axes_v[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) 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() bin_centers = 0.5 * (radial_bins[1:] + radial_bins[:-1]) @@ -593,7 +628,9 @@ class Plotter(Aggregator, BaseProcessor): fitlabel=None, **kwargs ): - + """ + Plot an histogram (PDF, etc ...) + """ if not ax_los is None: name = name + "_" + ax_los @@ -707,6 +744,9 @@ class Plotter(Aggregator, BaseProcessor): subname_y=None, **kwargs ): + """ + Generic plot routine, with name_x and name_y two path in the hdf5 file + """ if not node_arg is None: name_x, name_y = name_x + "_" + node_arg, name_y + "_" + node_arg @@ -856,12 +896,18 @@ class Plotter(Aggregator, BaseProcessor): hdf5_y.close() def _pspec(self, name, **kwargs): + """ + Plot power spectrum + """ del kwargs["run"] file_pspec = self.save.get_node("/hdf5/pspec").read() num = self.save.root._v_attrs.num getattr(pspec_read, "pspec_" + name)(file_pspec, ".", num, **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 yerr is None: (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) def overlay_kennicutt(self, n0, step): + """ + Add an overlay : kennicutt mass accretion + """ P.grid(False) ylim = P.ylim() (tmin, tmax) = P.xlim() @@ -934,6 +983,9 @@ class Plotter(Aggregator, BaseProcessor): P.ylim(ylim) def def_rules(self): + """ + That is where rules are defined + """ self.rules = { # Generic rules "plot": PlotRule( diff --git a/postprocessor.py b/postprocessor.py index 1b556cb..2d7879b 100644 --- a/postprocessor.py +++ b/postprocessor.py @@ -8,7 +8,7 @@ vol_func = lambda dset: dset["dx"] ** 3 # Volume function 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)) v_norm = np.sqrt(np.sum(dset["vel"] ** 2, axis=1)) # 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) -def _getter_B_int(dset): +def getter_B_int(dset): B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1)) return B_norm -def _getter_rho(dset): +def getter_rho(dset): return dset["rho"] @@ -302,27 +302,44 @@ class PostProcessor(HDF5Container): centers = 0.5 * (edges[1:] + edges[:-1]) 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() - B = _getter_B_int(self.cells) - rho = _getter_rho(self.cells) - B = np.log10(B) - rho = np.log10(rho) - abscisse = np.linspace(np.min(rho), np.max(rho), bins) + B = getter_B_int(self.cells) + rho = getter_rho(self.cells) + + if logbins: + 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) + + # For each cell, bin_number contains the number of the bins it belongs to 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) - 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)): B_mean[i] = np.mean(B[bin_number == i]) - values, edges = np.histogram(B, abscisse, weights=weights) - centers = 0.5 * (edges[1:] + edges[:-1]) - print("centers") + + # Get the center of each bin + 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: 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"]): 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) def _B_int(self, ax_los, z=0.0): + """ + Slice ont the intensity of the magnétic field + """ B_op = ScalarOperator( lambda dset: np.sqrt(np.sum(dset["Br"] ** 2, axis=1)), self._ro.info["unit_mag"], @@ -606,38 +626,6 @@ class PostProcessor(HDF5Container): pdf.attrs.var = np.var 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): csv_name = ( self.path @@ -829,17 +817,17 @@ class PostProcessor(HDF5Container): ), "cos_pdf": Rule( self, - partial(self._vol_pdf, _getter_abs_cos_vB, logbins=False), + partial(self._vol_pdf, getter_abs_cos_vB, logbins=False), "Global cos-PDF", "/hist", unit=cst.none, ), "Brho": Rule( self, - partial(self._Brho), - "Brho average", + self._Brho, + "Average of B as a function of rho", "/datasets", - unit=self.info["unit_mag"], + unit={"rho": self.info["unit_density"], "B": self.info["unit_mag"]}, ), # Profiles "axis": Rule(