From 8dfdf918a0aa361cd83266efcd54420e38e83266 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Tue, 21 Apr 2020 10:39:24 +0200 Subject: [PATCH] correct bugs --- comparator.py | 83 ++++++++++++++++++++++++++---------------------- plotter.py | 57 ++++++++++++++++++++++++++++----- postprocessor.py | 12 +++---- 3 files changed, 100 insertions(+), 52 deletions(-) diff --git a/comparator.py b/comparator.py index 4e2e5a2..9a5d128 100644 --- a/comparator.py +++ b/comparator.py @@ -137,39 +137,47 @@ class Comparator(Aggregator, HDF5Container): return np.array(prop) def _time_avg(self, name, start=None, end=None, span=None, group="/series"): - mean = np.zeros(len(self.runs)) - median = np.zeros(len(self.runs)) - std = np.zeros(len(self.runs)) - v_min = np.zeros(len(self.runs)) - v_max = np.zeros(len(self.runs)) - q975 = np.zeros(len(self.runs)) - q025 = np.zeros(len(self.runs)) - q16 = np.zeros(len(self.runs)) - q84 = np.zeros(len(self.runs)) + serie0 = self.save.get_node(group + "/" + name + "/" + self.runs[0]).read() + if len(serie0.shape) > 1: + shape = [len(self.runs)] + list(serie0.shape[1:]) + else: + shape = len(self.runs) + mean = np.zeros(shape) + median = np.zeros(shape) + std = np.zeros(shape) + v_min = np.zeros(shape) + v_max = np.zeros(shape) + q975 = np.zeros(shape) + q025 = np.zeros(shape) + q16 = np.zeros(shape) + q84 = np.zeros(shape) for i, run in enumerate(self.runs): serie = self.save.get_node(group + "/" + name + "/" + run).read() time = self.save.get_node(group + "/time/" + run).read() - mask = abs(serie) != np.inf + if len(serie.shape) <= 1: + mask = abs(serie) != np.inf - if not ((start, end, span) == (None, None, None)): - start_r, end_r = start, end - # Be sure that start_r and end_r are defined - if end_r is None: - if span is None or start_r is None: - end_r = time[-1] - else: - end_r = start_r + span - if start_r is None: - if span is None: - start_r = time[0] - else: - start_r = end_r - span + if not ((start, end, span) == (None, None, None)): + start_r, end_r = start, end + # Be sure that start_r and end_r are defined + if end_r is None: + if span is None or start_r is None: + end_r = time[-1] + else: + end_r = start_r + span + if start_r is None: + if span is None: + start_r = time[0] + else: + start_r = end_r - span + mask = ( + mask & (time >= start_r) & (time <= end_r) & np.isfinite(serie) + ) + serie = serie[mask] - mask = mask & (time >= start_r) & (time <= end_r) & np.isfinite(serie) - - mean[i] = np.mean(serie[mask], axis=0) - std[i] = np.std(serie[mask], axis=0) + mean[i] = np.mean(serie, axis=0) + std[i] = np.std(serie, axis=0) ( v_min[i], q025[i], @@ -178,7 +186,7 @@ class Comparator(Aggregator, HDF5Container): q84[i], q975[i], v_max[i], - ) = np.percentile(serie[mask], [0, 2.5, 16, 50, 84, 97.5, 100], axis=0) + ) = np.percentile(serie, [0, 2.5, 16, 50, 84, 97.5, 100], axis=0) return { "runs": self.runs, "mean": mean, @@ -403,16 +411,6 @@ class Comparator(Aggregator, HDF5Container): "q16": "16 percentile of {}".format(descr), "q84": "84 percentile of {}".format(descr), } - - # units={ - # "runs": cst.none, - # "mean": unit, - # "std": unit, - # "median": unit, - # "max": unit, - # "min": unit, - # "q025": unit, - # "q975": unit} units = unit if name is None: @@ -528,6 +526,14 @@ class Comparator(Aggregator, HDF5Container): group="/series", dependencies={"time": None, "rho_prof": "__parent__"}, ), + "time_coldens_pdf": Rule( + self, + partial( + self._time_series, partial(self.get_pp_value, "/hist/pdf_coldens") + ), + group="/series", + dependencies={"time": None, "pdf_coldens": "__parent__"}, + ), "time_pdf_slope_coldens": Rule( self, partial( @@ -560,6 +566,7 @@ class Comparator(Aggregator, HDF5Container): "time_pdf_slope_coldens", "turb_power", "time_rho_prof", + "time_coldens_pdf", ]: self._gen_rule_avg(name) diff --git a/plotter.py b/plotter.py index cb0c10d..042f6b1 100644 --- a/plotter.py +++ b/plotter.py @@ -126,7 +126,7 @@ class Plotter(Aggregator, BaseProcessor): name_full = name if rule.is_valid(arg): - if rule.kind == "classic": + if rule.kind == "classic" or rule.kind == "runs": try: runs = kwargs.pop("runs") if isinstance(runs, RunSelector): @@ -136,9 +136,17 @@ class Plotter(Aggregator, BaseProcessor): i = 0 for run in runs: files = [] - for num in self.nums[run]: + if rule.kind == "classic": + nums = self.nums[run] + else: + nums = [None] + for num in nums: plot_filename = self._find_filename(name_full, run, num) - save = tables.open_file(self.pp[run][num].filename) + + if rule.kind == "classic": + save = tables.open_file(self.pp[run][num].filename) + else: + save = tables.open_file(self.comp.filename, "r") try: self._plot_rule( rule, @@ -265,7 +273,10 @@ class Plotter(Aggregator, BaseProcessor): return r"{} = {}".format(prop_label, prop_value_str) if nml_key is None and label is None: - if "attrs" in self.save.root._v_attrs: + if ( + "attrs" in self.save.root._v_attrs + and run in self.save.root._v_attrs.attrs + ): label_run = r"{}".format(self.save.root._v_attrs.attrs[run].label) else: label_run = run @@ -479,8 +490,9 @@ class Plotter(Aggregator, BaseProcessor): def _plot_hist( self, name, - ax_los=None, - run="", + ax_los, + run, + group="/hist/", label=None, unit=None, unit_coeff=1.0, @@ -495,13 +507,15 @@ class Plotter(Aggregator, BaseProcessor): color=None, colors=None, nml_color=None, + fit=None, + fitlabel=None, **kwargs ): if not ax_los is None: name = name + "_" + ax_los - node = self.save.get_node("/hist/" + name) + node = self.save.get_node(group + name) if xlog is None: try: @@ -511,7 +525,12 @@ class Plotter(Aggregator, BaseProcessor): label, unit_old, unit = self._ax_label_unit(node, label, unit, unit_coeff) - values, centers = node.read() + if "mean" in node: + index = node["runs"].read().index(run) + values, centers = node["mean"].read()[index] + else: + values, centers = node.read() + if xlog: centers = centers + np.log10(unit_old.express(unit)) else: @@ -570,6 +589,11 @@ class Plotter(Aggregator, BaseProcessor): P.ylim([None, 1.0]) + if not fit is None: + self._overlay_fit( + centers, values, kind=fit, ls="--", lw=1.5, label=fitlabel + ) + def _plot( self, name_x, @@ -582,6 +606,7 @@ class Plotter(Aggregator, BaseProcessor): yunit=None, xunit_coeff=1.0, yunit_coeff=1.0, + ylog=False, fit=None, fitlabel=None, smooth=0, @@ -618,6 +643,9 @@ class Plotter(Aggregator, BaseProcessor): if grid: P.grid() + if ylog: + P.yscale("log") + if put_time: time = self.save.root._v_attrs.time * self.comp.info["unit_time"] time_str = self.pp_params.plot.time_fmt.format( @@ -878,6 +906,19 @@ class Plotter(Aggregator, BaseProcessor): "$\rho$-PDF", dependencies=["rho_pdf"], ), + "avg_coldens_pdf": PlotRule( + self, + partial( + self._plot_hist, + "avg_time_coldens_pdf_z", + group="/comp/", + xlog=True, + put_time=False, + ), + "PDF", + kind="runs", + dependencies={"avg_time_coldens_pdf": "z"}, + ), "T_pdf": PlotRule( self, partial(self._plot_hist, "T_pdf"), "T-PDF", dependencies=["T_pdf"] ), diff --git a/postprocessor.py b/postprocessor.py index a9fa573..8323519 100644 --- a/postprocessor.py +++ b/postprocessor.py @@ -75,7 +75,7 @@ class PostProcessor(HDF5Container): # Set the extend of the image self._radius = 0.5 / self.pp_params.pymses.zoom - self._lbox = self.info["boxlen"] + self.lbox = self.info["boxlen"] center = self.pp_params.pymses.center im_extent = [ (-self._radius + center[0]), @@ -91,7 +91,7 @@ class PostProcessor(HDF5Container): self.save.root._v_attrs.dir = os.path.dirname(path) self.save.root._v_attrs.run = os.path.basename(path) self.save.root._v_attrs.num = num - self.save.root._v_attrs.lbox = self._lbox + self.save.root._v_attrs.lbox = self.lbox self.save.root._v_attrs.unit_length = self.info["unit_length"] self.save.root._v_attrs.time = time @@ -361,7 +361,7 @@ class PostProcessor(HDF5Container): # Operator to compute the angular speed times rho def omega_rho_func(dset): pos = dset.get_cell_centers() - pos = pos - (self.pp_params.disk.pos_star / self._lbox) + pos = pos - (self.pp_params.disk.pos_star / self.lbox) xx = pos[:, :, 0] yy = pos[:, :, 1] rc = np.sqrt(xx ** 2 + yy ** 2) # cylindrical radius @@ -401,7 +401,7 @@ class PostProcessor(HDF5Container): dmap_cs = rt_cs.process(self._cam[ax_los]) dmap_col = self.save.root.maps.coldens_z.read() map_Q = ( - (self._lbox * dmap_cs.map.T) + (self.lbox * dmap_cs.map.T) * dmap_omega.map.T / (np.pi * self.G * dmap_col) ) @@ -412,7 +412,7 @@ class PostProcessor(HDF5Container): Computes radial bins (for disk) """ pos_star = self.pp_params.disk.pos_star - im_extent = self.save.root.maps._v_attrs.im_extent + im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * self.lbox # radius of the corner of the box plus a margin rad_of_box = ( @@ -442,7 +442,7 @@ class PostProcessor(HDF5Container): """ Computes the radius from the center """ - im_extent = self.save.root.maps._v_attrs.im_extent + im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * self.lbox map_size = self.pp_params.pymses.map_size pos_star = self.pp_params.disk.pos_star