From ed55a0cba19ef00798dae4ee1404026a994bfe75 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Mon, 4 Nov 2019 16:05:09 +0100 Subject: [PATCH] Began the work on a new pipeline using the new HDF5 version --- pipeline.py | 284 +++++++++++++++++ plotter.py | 353 +++++++++++++++------ postprocessor.py | 642 +++++++++++++++++++++++++++++++------- params.py => pp_params.py | 24 +- 4 files changed, 1083 insertions(+), 220 deletions(-) create mode 100644 pipeline.py rename params.py => pp_params.py (56%) diff --git a/pipeline.py b/pipeline.py new file mode 100644 index 0000000..3f0a2ae --- /dev/null +++ b/pipeline.py @@ -0,0 +1,284 @@ +# coding: utf-8 + +import os +import glob +from shutil import copy +import argparse +import time +import numpy as np +from functools import reduce + +from pp_params import * +from plotter import * +from postprocessor import * + +fake_pp = PostProcessor() + +parser = argparse.ArgumentParser() + +input_args = parser.add_argument_group('input', "Input selection") +input_args.add_argument("runs", help='name of runs', nargs='*', + default=[]) + +input_args.add_argument("-ip", "--input_path", + help="specify input directory", + default="/home/nbrucy/simus/") +input_args.add_argument("-p", "--project", + help="specify project name (directory within the input directory)", + default="disk") + + +input_args.add_argument("-wo", "--which_inputs", + choices=['all', 'id', 'time'], + help="Select inputs by time range, id range or all of them", + default='all') +input_args.add_argument("-b", "--begin", + help="id of first input", + type=int, + default=1) +input_args.add_argument("-e", "--end", + help="id of last input", + type=int, + default=100) +input_args.add_argument("-s", "--step", + help="step between two input", + type=int, + default=1) +input_args.add_argument("-tb", "--time_begin", + help="time of first input", + type=float, + default=0.) +input_args.add_argument("-te", "--time_end", + help="time of last input", + type=float, + default=6.) + +input_args.add_argument("-w", "--watch", + help="wait and watch for missing inputs", + action='store_true') +input_args.add_argument("--skip", + help="skip failed loadings", + action='store_true') +input_args.add_argument("-wt", "--waiting_time", + help="time between to successive try when watching new inputs (in second)", + type=int, + default=120) +input_args.add_argument("-af", "--allowed_failures", + help="number of allowed failures when waiting", + default=30) + +output_args = parser.add_argument_group('output', "Output configuration") + +output_args.add_argument("-op", "--output_path", + help="specify output directory", + default='/home/nbrucy/visus/') + +output_args.add_argument("--tag", + help="Add a special tag on output filemanes", + default='') + +output_args.add_argument("-owr", "--overwrite", + help="Overwrite outputs", + action='store_true') + + +output_args.add_argument("-owrd", "--overwrite_dependencies", + help="Overwrite outputs for dependencies", + action='store_true', + default=None) + +pp_args = parser.add_argument_group('postproc', "Post Processing configuration") + + +pp_args.add_argument("--process", + help="Individual rules to apply", + choices=fake_pp.rules.keys(), + default=[], + nargs='*') + +pp_args.add_argument("-pargs", "--process_args", + help="Args to give to process rules", + default=['x', 'y', 'z'], + nargs='*') + +pp_args.add_argument("--compare", + help="Time and inter run comparaison", + default=[], + nargs='*') + +pp_args.add_argument("-cargs", "--compare_args", + help="Args to give to process rules", + default=[None], + nargs='*') + +pp_args.add_argument("--plot", + help="Plot rules", + default=[], + nargs='*') + +pp_args.add_argument("-plargs", "--plot_args", + help="Args to give to plot rules", + default=['x', 'y', 'z'], + nargs='*') + +pp_args.add_argument("-d", "--disk", + help="Specify this for disk simulations", + action='store_true') +pp_args.add_argument("--fft", + help="use quick and dirty fft rendering", + action='store_true') +pp_args.add_argument("--zoom", + help="zoom", + type=float, + default=2.) +pp_args.add_argument("-ms", "--map_size", + help="size of the maps created in he map mode (in pixel)", + type=int, + default=1024) + +pp_args.add_argument("--nb_bin", + help="Number of bins for azimuthal averages", + type=int, + default=50) +pp_args.add_argument("--pdf_nb_bin", + help="Number of bins for PDF", + type=int, + default=50) +pp_args.add_argument("--binning", + help="Kind of binning (logarithmic or linear)", + choices=['log', 'lin'], + default='log') + +plot_args = parser.add_argument_group('plot', "Plot configuration") + +plot_args.add_argument("--colormap", + help="Colormap used", + choices=P.colormaps(), + default='plasma') +plot_args.add_argument("--format", + help="Format of the plot images", + choices=['png', 'jpeg', 'pdf', 'svg', 'ps'], + default='jpeg') +plot_args.add_argument("--dpi", + help="Resolution of the plot images", + type=int, + default=400) +plot_args.add_argument("--beamer", + help="Beamer mode", + action='store_true') + +args = parser.parse_args() + +project = args.project +runs = args.runs +storage_in = args.input_path +storage_out = args.output_path + +pp_params = Params() + +pp_params.out.zoom = args.zoom +pp_params.out.tag = args.tag +pp_params.out.map_size = args.map_size + +pp_params.pymses.fft = args.fft + + +pp_params.disk.on = args.disk +pp_params.disk.binning = args.binning +pp_params.disk.nb_bin = args.nb_bin + +pp_params.pdf.nb_bin = args.pdf_nb_bin + +# extension for out files +P.style.use("seaborn-deep") +if args.format == 'pdf': + P.style.use("~/.config/matplotlib/pdf.mplstyle") + +if args.beamer: + P.rcParams['font.family'] = 'sans-serif' + P.rcParams['figure.figsize'] = (7, 4.5) + +# Plot properties +P.rcParams['image.cmap'] = args.colormap +P.rcParams['savefig.dpi'] = args.dpi +P.rcParams['lines.linewidth'] = 2 +P.rcParams['lines.markersize'] = 10 +P.rcParams["errorbar.capsize"] = 4 + +# List of id that were successfully computed +nums_success = dict() + + +# Go through all runs +for run in runs: + path_suffix = project + '/' + run + path_in = storage_in + path_suffix + path_out = storage_out + path_suffix + + if args.tag == '': + tag_run = run + else: + tag_run = run + '_' + args.tag + + if not os.path.exists(path_out): + os.makedirs(path_out) + try: + copy(path_in + '/disk.nml', path_out) + copy(path_in + '/output_00001/compilation.txt', path_out) + except: + pass + + nums_success[run] = [] + + if args.which_inputs in ['all', 'time'] : + names = glob.glob(path_in + "/output_[0-9][0-9][0-9][0-9][0-9]") + nums_all = [int(n.split('/')[-1].split('_')[1]) for n in names] + nums_all = np.sort(nums_all) + if args.which_inputs == 'all': + nums = nums_all + else: + time = [get_time(path_in, n) for n in nums_all] + nums = [n for i,n in enumerate(nums_all) if time[i] >= args.time_begin + and time[i] < args.time_end] + else: + nums = range(args.begin, args.end + 1, args.step) + + for num in nums: + failures = 0 + success = False + + while not success: + try: + if len(args.process) > 0 and len(args.plot) > 0: + pp = PostProcessor(run, num, pp_params=pp_params) + pp.process(args.process, args.process_args, + overwrite=args.overwrite, overwrite_dep=args.overwrite_dependencies) + + pltter = Plotter(filename=pp.filename) + pltter.plot(args.plot, args.plot_args, + overwrite=args.overwrite, overwrite_dep=args.overwrite_dependencies) + + # If we are here, success ! + success = True + nums_success[run].append(num) + except (ValueError, IOError, KeyError) as e: + print(e) + if(args.watch and failures < args.allowed_failures): + failures = failures + 1 + print("ERROR: Unable to proceed for run {} output {}.\ + Trying again in {} s ({} tries remaining)".format(run, num, + args.waiting_time, args.allowed_failures - failures)) + time.sleep(args.waiting_time) + elif args.skip: + break + else: + raise + +if len(args.compare) > 0: + path_in = storage_in + project + path_out = storage_out + project + + cc = Comparator(path_in, runs, nums_success, path_out=path_out, pp_params=pp_params) + cc.process(args.compare, args.compare_args, + overwrite=args.overwrite, overwrite_dep=args.overwrite_dependencies) + diff --git a/plotter.py b/plotter.py index 46dd738..091a97b 100644 --- a/plotter.py +++ b/plotter.py @@ -9,11 +9,15 @@ if os.environ.get('DISPLAY','') == '': print('No display found. Using non-interactive Agg backend') mpl.use('Agg') from matplotlib.patches import Polygon -import pylab as Pfrom scipy.stats import linregress +import pylab as P +from scipy.stats import linregress import matplotlib.patches as mpatches from matplotlib.collections import PatchCollection +from functools import partial +from numpy.polynomial.polynomial import polyfit from pp_params import * +from postprocessor import Rule, BaseProcessor P.rcParams['image.cmap']='plasma' @@ -23,7 +27,22 @@ tex_params= {'text.latex.preamble' : [r'\usepackage{amsmath}']} P.rcParams.update(tex_params) -class Plotter: +class PlotRule(Rule): + + def plot(self, arg): + return self.process_fn(arg) + + def is_valid(self, arg): + save = self.postproc.save + valid = True + for dep in self.dependencies: + if not arg is None: + valid = valid and dep + '_' + str(arg) in save + else: + valid = valid and dep in save + return arg in self.args_ok and valid and self.is_valid_add(arg) + +class Plotter(BaseProcessor): """ This class loads derived quantities and plot them """ @@ -36,7 +55,7 @@ class Plotter: G = 1. # Gravitational constant - def __init__(self, path_out='.', filename=None, pp_params=Params()): + def __init__(self, filename=None, path_out='.', num=None, pp_params=Params()): self.pp_params = pp_params @@ -60,132 +79,264 @@ class Plotter: else: self.filename = filename - def plot_list(self, to_plot_list, axes, overwrite=False): - self._file_out = tables.open_file(self.filename, mode="r") - maps = self._file_out.root.maps + self.log_id = "[plot {}] ".format(self.pp_params.out.tag) + self.def_rules() - for ax_los in axes: - for name in to_plot_list: - name_full = name + '_' + ax_los - plot_filename = self._find_filename(name_full) - if overwrite or not os.path.exists(plot_filename): - self._plot_map(name, ax_los) - P.savefig(plot_filename) - else: - print("Data for {} is already computed, skipping...".format(name_full)) + def plot(self, to_plot_list, args=[None], overwrite=False): + """ + Plot the data in to_plot_list and save them + """ + self.process(to_plot_list, args, overwrite, False) - self._file_out.close() + def _process_single(self, name, rule, arg, overwrite=False, overwrite_dep=False, just_done=[]): + done = self._plot_rule(name, rule, arg, overwrite) + return [] + + def _plot_rule(self, name, rule, arg, overwrite): + if not arg is None: + name_full = rule.group + '/' + name + '_' + str(arg) + else: + name_full = rule.group + '/' + name + + if rule.is_valid(arg): + plot_filename = self._find_filename(name_full) + if overwrite or not os.path.exists(plot_filename): + rule.plot(arg) + P.tight_layout(pad=1) + P.savefig(plot_filename) + P.close() + self._log("{} plotted".format(name_full), "SUCCESS") + else: + self._log("Plot {} is already done, skipping...".format(name_full)) + else: + self._log("{} is not valid in this context".format(name_full), "ERROR") + + + def _find_filename(self, name_full): + if not self.pp_params.out.tag == '': + tag_name = '_' + self.pp_params.out.tag + else : + tag_name = '' + + if 'num' in self.save.root._v_attrs: + num = self.save.root._v_attrs.num + return (self.path_out + '/' + name_full + + tag_name + '_' + format(num,'05') + + self.pp_params.plot.out_ext) + else: + return self.path_out + '/' + name_full + tag_name + self.pp_params.plot.out_ext - def _plot_map(self, name, ax_los): + def _plot_map(self, name, ax_los, label=None, cmap='plasma', vmin=None, vmax=None, overlays=[]): P.figure() - ax_h = self._axes_h[ax_los] ax_v = self._axes_v[ax_los] - im_extent = self._file_out.root.maps._v_attrs.im_extent - radius = self._file_out.root.maps._v_attrs.radius - center = self._file_out.root.maps._v_attrs.center + im_extent = self.save.root.maps._v_attrs.im_extent - if (name == 'Q' and not ax_los == 'z') or name == 'levels' or name=='speed': - return + dmap = self.save.get_node('/maps/{}_{}'.format(name, ax_los)).read() - dmap = self._file_out.get_node('/maps/{}_{}'.format(name, ax_los)).read() - - if name == 'Q' : - im = P.imshow(dmap, - extent=im_extent, - origin='lower', - cmap='RdBu', - norm=mpl.colors.LogNorm(), - vmin=0.01, vmax=100.) - elif name == 'jeans_ratio' : - im = P.imshow(dmap, - extent=im_extent, - origin='lower', - cmap='RdBu', - norm=mpl.colors.LogNorm(), - vmin=0.1, vmax=1000.) - else: - im = P.imshow(dmap, - extent=im_extent, - origin='lower', - norm=mpl.colors.LogNorm()) - - P.locator_params(axis=ax_h, nbins=pp.params.plot.ntick) - P.locator_params(axis=ax_v, nbins=pp.params.plot.ntick) - - if(self.pp_params.put_title): - pass + im = P.imshow(dmap, + extent=im_extent, + origin='lower', + cmap=cmap, + norm=mpl.colors.LogNorm()) + im.set_clim(vmin, vmax) + P.locator_params(axis=ax_h, nbins=self.pp_params.plot.ntick) + P.locator_params(axis=ax_v, nbins=self.pp_params.plot.ntick) P.xlabel(self._ax_title[ax_h]) P.ylabel(self._ax_title[ax_v]) cbar = P.colorbar(im) - if name == 'coldens': - cbar.set_label(r'$\Sigma$ (code)') + if not label is None: + cbar.set_label(label) - if pp.params.set_lim: - im.set_clim(0.01, 100) + for plot_overlay in overlays: + plot_overlay(ax_los) - # if 'levels' in names: - # map_level = self._file_out.get_node('/maps/{}_{}'.format('levels', ax_los)).read() - # # Computing linewidths - # levels_ar = np.arange(np.min(map_level), np.max(map_level) + 1) - # lw = np.ones(levels_ar.size) * 2 - # lvl_th = 8 # Level threeshold for reducing linewidths - # lw[levels_ar >= lvl_th] = lw[levels_ar >= lvl_th]**(lvl_th - levels_ar[levels_ar >= lvl_th]) - # lw[levels_ar < lvl_th] = 1. + def _overlay_levels(self, ax_los): + 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) + lw = np.ones(levels_ar.size) * 2 + lvl_th = 8 # Level threeshold for reducing linewidths + lw[levels_ar >= lvl_th] = lw[levels_ar >= lvl_th]**(lvl_th - levels_ar[levels_ar >= lvl_th]) + lw[levels_ar < lvl_th] = 1. - # cont = P.contour(map_level, - # extent=im_extent, - # origin='lower', - # colors='white', - # linewidths=lw, - # levels=levels_ar) - # cont.levels = cont.levels + 1 + cont = P.contour(map_level, + extent=self.save.root.maps._v_attrs.im_extent, + origin='lower', + colors='grey', + linewidths=lw, + levels=levels_ar) + cont.levels = cont.levels + 1 - # P.clabel(cont, - # cont.levels[cont.levels < 11], - # inline=1, fontsize=8., fmt='%1d') - elif name == 'rho': - cbar.set_label(r'$\rho$ (code)') + P.clabel(cont, + cont.levels[cont.levels < 11], + inline=1, fontsize=8., fmt='%1d'); - if 'speed' in names: - dmap_vh = self._file_out.get_node('/maps/{}{}_{}'.format('v', ax_h, ax_los)).read() - dmap_vv = self._file_out.get_node('/maps/{}{}_{}'.format('v', ax_v, ax_los)).read() + def _overlay_speed(self, ax_los): + ax_h = self._axes_h[ax_los] + ax_v = self._axes_v[ax_los] + dmap_vh = self.save.get_node('/maps/speed_h_{}'.format(ax_los)).read() + dmap_vv = self.save.get_node('/maps/speed_v_{}'.format(ax_los)).read() - vel_red = self.pp_params.plot.vel_red + vel_red = self.pp_params.plot.vel_red + radius = self.save.root.maps._v_attrs.radius + center = self.save.root.maps._v_attrs.center + lbox = self.save.root._v_attrs.lbox - map_vh_red = dmap_vh[::vel_red,::vel_red] # take only a subset of velocities - map_vv_red = dmap_vv[::vel_red,::vel_red] - nh = map_vh_red.shape[0] - nv = map_vv_red.shape[1] - vec_h = (np.arange(nh)*2./nh*radius - radius + center[0] + radius/nh) * lbox - vec_v = (np.arange(nv)*2./nv*radius - radius + center[1] + radius/nv) * lbox - hh, vv = np.meshgrid(vec_h,vec_v) - max_v = np.max(np.sqrt(map_vh_red**2 + map_vv_red**2)) + map_vh_red = dmap_vh[::vel_red,::vel_red] # take only a subset of velocities + map_vv_red = dmap_vv[::vel_red,::vel_red] + nh = map_vh_red.shape[0] + nv = map_vv_red.shape[1] + vec_h = (np.arange(nh)*2./nh*radius - radius + center[0] + radius/nh) * lbox + vec_v = (np.arange(nv)*2./nv*radius - radius + center[1] + radius/nv) * lbox + hh, vv = np.meshgrid(vec_h,vec_v) + max_v = np.max(np.sqrt(map_vh_red**2 + map_vv_red**2)) - Q = P.quiver(hh, vv, map_vh_red, map_vv_red, units='width') - P.quiverkey(Q, 0.7, 0.95, max_v, - r'$'+str(max_v)[0:4]+'$ (code)', labelpos='E', coordinates='figure') + Q = P.quiver(hh, vv, map_vh_red, map_vv_red, units='width', color='grey') + P.quiverkey(Q, 0.6, 0.98, max_v, + r'$'+str(max_v)[0:4]+'$ (code)', labelpos='E', coordinates='figure') - elif name == 'T': - cbar.set_label(r'$T (code)$') - elif name == 'Q': - cbar.set_label(r'$Q$') - elif name == 'jeans': - cbar.set_label(r'Jeans\'s lenght') + def _plot_radial(self, name, ax_los='z', label=None, xlog=False, ylog=False): + P.figure() + radial_bins = self.save.get_node('/radial/radial_bins_' + ax_los).read() + bin_centers = 0.5 * (radial_bins[1:] + radial_bins[:-1]) + mean_bin = self.save.get_node('/radial/{}_{}'.format(name, ax_los)).read() + + P.grid() + P.xlabel(r'$r$') + + if xlog: + P.xscale('log') + if ylog: + P.yscale('log') + P.plot(bin_centers, mean_bin) + + if not label is None: + P.ylabel(label) + + def _plot_hist(self, name, ax_los='z', label=None, ylog=False): + P.figure() + pdf = self.save.get_node('/hist/' + name + '_' + ax_los) + values, centers = pdf.read() + width = centers[1] - centers[0] + P.bar(centers, values, width, log=ylog) + P.grid() + if not label is None: + P.xlabel(label) + + if '/hist/fit_' + name + '_' + ax_los in self.save: + slope = pdf.attrs.slope + origin = pdf.attrs.origin + P.plot(centers, 10**(slope*centers + origin), '--', linewidth=2, color='orange') + + P.ylim([None, 1.]) + + def _plot(self, name_x, name_y, xlabel=None, ylabel=None, linearfit=False): + P.figure() + node_x = self.save.get_node(name_x) + node_y = self.save.get_node(name_y) + + if xlabel is None: + if 'label' in node_x: + xlabel = name_x._vattrs.label + else: + xlabel = name_x + + if ylabel is None: + if 'label' in node_y: + ylabel = name_y._vattrs.label + else: + ylabel = name_y + + x = node_x.read() + if node_y._v_attrs.CLASS == 'ARRAY': + y = node_y.read() + P.plot(x, y, fmt='*') else: - cbar.set_label(name) + y = node_y.mean.read() + yerr = node_y.std.read() + P.errorbar(x, y, yerr=y, fmt='*') - def _find_filename(self, name_full): - num = self._file_out.root._v_attrs.num - return (self.path_out + '/' + name_full + '_' + - self.pp_params.out.tag + '_' + format(num,'05') + - self.pp_params.plot.out_ext) + P.xlabel(xlabel) + P.ylabel(ylabel) + if linearfit: + if node_y._v_attrs.CLASS == 'ARRAY': + (a, b, rho, _, stderr) = linregress(node_x.read(), node_y.read()) + else: + c = polyfit(node_x.read(), node_y.mean.read(), 1, + w = [1.0 / ty for ty in node_y.std.read()], full=False) + b, a = c[0], c[1] + + P.plot(x, a*y + b, '--', linewidth=1.5) + + P.grid() + + + + def def_rules(self): + self.rules = { + 'coldens' : PlotRule(self, lambda ax_los: + self._plot_map('coldens', ax_los, label=r'$\Sigma$ (code)', vmin=0.01, vmax=100), + "Column density", ['/maps/coldens']), + 'rho' : PlotRule(self, lambda ax_los: + self._plot_map('rho', ax_los, label=r'$\rho$ (code)'), + "Density slice", ['/maps/rho']), + 'coldens_l' : PlotRule(self, lambda ax_los: + self._plot_map('coldens', ax_los, label=r'$\Sigma$ (code)', + vmin=0.01, vmax=100, overlays=[self._overlay_levels]), + "Column density", ['/maps/coldens', '/maps/levels']), + 'rho_v' : PlotRule(self, lambda ax_los: + self._plot_map('rho', ax_los, label=r'$\rho$ (code)', + overlays=[self._overlay_speed]), + "Density slice", ['/maps/rho', '/maps/speed_h', '/maps/speed_v']), + 'jeans_ratio' : PlotRule(self, lambda ax_los: + self._plot_map('jeans_ratio', ax_los, vmin=0.1, vmax=100, + cmap='RdBu_r', + overlays=[self._overlay_levels]), + "Jeans' lenght divided by the max resolution", + dependencies=['/maps/jeans_ratio', '/maps/levels']), + 'Q' : PlotRule(self, lambda ax_los: + self._plot_map('rho', ax_los, label=r'$Q$', vmin=0.01, vmax=100, cmap='RdBu_r'), + "Toomre Q parameter for a Keplerian disk", + dependencies=['/maps/Q'], args_ok=['z']) + + } + + averageables = ['coldens', 'rho', 'T', 'Q'] + for name in averageables: + self.rules['rad_' + name] = PlotRule(self, partial(self._plot_radial, 'rad_avg_' + name, + label=name, xlog=True, ylog=True), + "Azimuthal average of {}".format(name), + dependencies=['/radial/radial_bins', '/radial/rad_avg_' + name], + args_ok=['z']) + + self.rules['fluct_' + name] = PlotRule(self, partial(self._plot_map, 'fluct_' + name, + vmin=0.01, vmax=100, cmap='RdBu_r', + label='{}/avg({})'.format(name, name)), + "Fluctuation wrt to average of {}".format(name), + dependencies=['/maps/fluct_' + name], + args_ok=['z']) + self.rules['pdf_' + name] = PlotRule(self, partial(self._plot_hist, 'pdf_' + name, ylog=True, + label='{}/avg({})'.format(name, name)), + "Probability density function of {} fluctuations".format(name), + dependencies=['/hist/pdf_' + name], + args_ok=['z']) + + + self.rules.update({ + 'kappa_beta' : PlotRule(self, partial(self._plot, '/comp/beta', '/comp/avg_pdf_slope_coldens', + linearfit=True), + args_ok=[None], dependencies=['/comp/beta', '/comp/avg_pdf_slope_coldens']), + 'sink_mass' : PlotRule(self, partial(self._plot, '/series/time', '/series/sinks_mass', + linearfit=True), + args_ok=[None], dependencies=['/series/time', '/series/sinks_mass']) + }) class InteractiveGUI: diff --git a/postprocessor.py b/postprocessor.py index c75e676..f9c47f5 100644 --- a/postprocessor.py +++ b/postprocessor.py @@ -12,30 +12,135 @@ from pymses.sources.hop.file_formats import * from pymses.analysis import Camera, raytracing, slicing, splatting from pymses.filters import CellsToPoints from pymses.analysis import ScalarOperator, FractionOperator, MaxLevelOperator +from scipy.stats import linregress +from functools import partial +from abc import ABCMeta, abstractmethod from pp_params import * class Rule: - - def __init__(self, process, description, group='', dependencies=[], axes=['x', 'y', 'z'], - is_valid=lambda save, ax:True): + def __init__(self, postproc, process, description='', group='', dependencies=[], args_ok=['x', 'y', 'z'], + is_valid=lambda arg:True): + self.postproc = postproc self.process_fn = process self.dependencies = dependencies self.is_valid_add = is_valid self.group = group - self.axes = axes + self.args_ok = args_ok self.description = description - def process(self, ax_los): - return self.process_fn(ax_los) + def process(self, arg): + if not arg is None: + return self.process_fn(arg) + else: + return self.process_fn() - def is_valid(self, save, ax): + def is_valid(self, arg): + save = self.postproc.save valid = True for dep in self.dependencies: - valid = valid and self.group + '/' + dep + '_' + ax in save - return ax in self.axes and valid and self.is_valid_add(save, ax) + rule_dep = self.postproc.rules[dep] + if not arg is None: + valid = valid and rule_dep.group + '/' + dep + '_' + str(arg) in save + else: + valid = valid and rule_dep.group + '/' + dep in save + return arg in self.args_ok and valid and self.is_valid_add(arg) -class PostProcessor: + +class BaseProcessor: + """ + Base class for processors, should not be instanciated + """ + __metaclass__ = ABCMeta + + @abstractmethod + def __init__(self): + self.def_rules() + + log_id = "" + + def _log(self, string, status=""): + if len(status) > 0: + print(status + ": " + self.log_id + string) + else: + print(self.log_id + string) + + def process(self, to_process_list, args=[None], overwrite=False, overwrite_dep=None): + """ + Render the data in to_process_list and save them + """ + if overwrite_dep is None: + overwrite_dep = overwrite + + self.overwrite_dep = overwrite_dep + just_done = [] # Computations done within this call + + self.save = tables.open_file(self.filename, mode="a") + for name in to_process_list: + if name in self.rules: + rule = self.rules[name] + for arg in args: + just_done = self._process_single(name, rule, arg, overwrite, just_done) + else: + self._log("{} is unknown, allowed rules are {}".format(name, self.rules.keys()), "ERROR") + self.save.close() + + def _process_single(self, name, rule, arg, overwrite=False, just_done=[]): + # Solve dependencies + for dep in rule.dependencies: + if dep in self.rules: + rule_dep = self.rules[dep] + just_done = self._process_single(dep, rule_dep, arg, self.overwrite_dep, just_done) + else: + self._log("Dependency {} for {} is unknown".format(dep, name), "ERROR") + # Process rule + done = self._process_rule(name, rule, arg, overwrite, just_done) + return just_done + [done] + + def _process_rule(self, name, rule, arg, overwrite=False, just_done=[]): + if not arg is None: + name_full = rule.group + '/' + name + '_' + str(arg) + else: + name_full = rule.group + '/' + name + + if rule.is_valid(arg): + if not name_full in just_done: + if overwrite or not name_full in self.save: + self._log("Processing {}".format(name_full)) + data = rule.process(arg) + self._save_data(name_full, data, rule.description) + self._log("Data for {} computed".format(name_full), "SUCCESS") + return name_full + else: + self._log("Data for {} is already computed, skipping...".format(name_full)) + else: + self._log("{} is not valid in this context".format(name_full), "ERROR") + + + def _save_data(self, name_full, data, description): + """ + Save data in the HDF5 structure, overwrite if necessary + """ + if name_full in self.save: + self.save.remove_node(name_full, recursive=True) + + if not type(data) == dict: + self.save.create_array(os.path.dirname(name_full), os.path.basename(name_full), + data, description, createparents=True) + else: + for key in data: + if type(description) == dict: + self.save.create_array(name_full, key, + data[key], description[key], createparents=True) + else: + self.save.create_array(name_full, key, + data[key], description, createparents=True) + + @abstractmethod + def def_rules(self): + pass + +class PostProcessor(BaseProcessor): """ This class enable to compute and save derived quantities from the raw output """ @@ -47,124 +152,104 @@ class PostProcessor: G = 1. # Gravitational constant - - def __init__(self, path, num, path_out=None, pp_params=Params()): + def __init__(self, path=None, num=None, path_out=None, filename=None, pp_params=Params()): """ Creates the basic structures needed for the outputs """ - # TODO : Make possible to load the HDF5 file even without the original file - self.pp_params = pp_params + if not path is None and not num is None: + # TODO : Make possible to load the HDF5 file even without the original file + self.pp_params = pp_params - # Determining output directory - if (path_out is None): - path_out = path + # Determining output directory + if (path_out is None): + path_out = path - # Open outfile - if not pp_params.out.tag == '': - tag_name = '_' + pp_params.out.tag - else : - tag_name = '' + # Open outfile + if not pp_params.out.tag == '': + tag_name = pp_params.out.tag + '_' + else : + tag_name = '' - self.filename = (path_out + '/postproc_' + - tag_name + format(num,'05') + '.h5') - self.save = tables.open_file(self.filename, mode="a", - title=os.path.basename(path) + format(num,'05')) + self.filename = (path_out + '/postproc_' + + tag_name + format(num,'05') + '.h5') + self.save = tables.open_file(self.filename, mode="a", + title=os.path.basename(path)+ '_' + format(num,'05')) - # Ramses Output - self._ro = pymses.RamsesOutput(path, num, order=pp_params.pymses.order) - self._amr = self._ro.amr_source(["rho","vel","P"]) + # Ramses Output + self.path = path + self.run = os.path.basename(path) + self.num = num + self._ro = pymses.RamsesOutput(path, num, order=pp_params.pymses.order) + self._amr = self._ro.amr_source(["rho","vel","P"]) - # Density operator - self._rho_op = ScalarOperator(lambda dset: dset["rho"], self._ro.info["unit_density"]) + # Density operator + self._rho_op = ScalarOperator(lambda dset: dset["rho"], self._ro.info["unit_density"]) - # Density ray tracer - if(pp_params.pymses.fft): - self._rt = splatting.SplatterProcessor(self._amr, self._ro.info, self._rho_op) - else: - self._rt = raytracing.RayTracer(self._amr, self._ro.info, self._rho_op) + # Density ray tracer + if(pp_params.pymses.fft): + self._rt = splatting.SplatterProcessor(self._amr, self._ro.info, self._rho_op) + else: + self._rt = raytracing.RayTracer(self._amr, self._ro.info, self._rho_op) - # Set the extend of the image - self._radius = 0.5 / pp_params.out.zoom - self._lbox = self._ro.info['boxlen'] - center = pp_params.out.center - im_extent = [(- self._radius + center[0]) * self._lbox, - ( self._radius + center[0]) * self._lbox, - (- self._radius + center[1]) * self._lbox, - ( self._radius + center[1]) * self._lbox] + # Set the extend of the image + self._radius = 0.5 / pp_params.out.zoom + self._lbox = self._ro.info['boxlen'] + center = pp_params.out.center + im_extent = [(- self._radius + center[0]) * self._lbox, + ( self._radius + center[0]) * self._lbox, + (- self._radius + center[1]) * self._lbox, + ( self._radius + center[1]) * self._lbox] - # Get time - time = self._ro.info['time'] # time in codeunits + # Get time + time = self._ro.info['time'] # time in codeunits - # Set post processing attributes - self.save.root._v_attrs.num = num - self.save.root._v_attrs.lbox = self._lbox - self.save.root._v_attrs.time = time + # Set post processing attributes + 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.time = time - self.save.root.maps._v_attrs.center = center - self.save.root.maps._v_attrs.radius = self._radius - self.save.root.maps._v_attrs.im_extent = im_extent + if not '/maps' in self.save: + self.save.create_group('/', 'maps', '2D maps') + self.save.root.maps._v_attrs.center = center + self.save.root.maps._v_attrs.radius = self._radius + self.save.root.maps._v_attrs.im_extent = im_extent - # Initialize cameras - self._cam = dict() - for ax_los in self._ax_nb : # los = line of sight - ax_h = self._axes_h[ax_los] - ax_v = self._axes_v[ax_los] + # Initialize cameras + self._cam = dict() + for ax_los in self._ax_nb : # los = line of sight + ax_h = self._axes_h[ax_los] + ax_v = self._axes_v[ax_los] + + self._cam[ax_los] = Camera(center=pp_params.out.center, + line_of_sight_axis=ax_los, + region_size=[2.*self._radius, 2.*self._radius], + distance=self._radius, + far_cut_depth=self._radius, + up_vector=ax_v, + map_max_size=pp_params.out.map_size) + + self._add_metadata() + self.save.close() + + self.log_id = "[{}, {}] ".format(self.run, self.num) - self._cam[ax_los] = Camera(center=pp_params.out.center, - line_of_sight_axis=ax_los, - region_size=[2.*self._radius, 2.*self._radius], - distance=self._radius, - far_cut_depth=self._radius, - up_vector=ax_v, - map_max_size=pp_params.out.map_size) - self.save.close() self.def_rules() - def process(self, to_process_list, axes, overwrite=False): + def _add_metadata(self): """ - Render the data in to_process_list and save them + Add additional metadata to the file """ - self.save = tables.open_file(self.filename, mode="a") - for name in to_process_list: - if name in self.rules: - rule = self.rules[name] - for ax_los in axes: - # Solve dependencies - for dep in rule.dependencies: - if dep in self.rules: - rule_dep = self.rules[dep] - self._process_rule(dep, rule_dep, ax_los, overwrite) - else: - print("ERROR: Dependency {} for {} is unknown".format(dep, name)) - # Process rule - self._process_rule(name, rule, ax_los, overwrite) + + # Beta for the beta cooling + if not (self.pp_params.disk.beta is None or self.pp_params.disk.beta == False): + if type(self.pp_params.disk.beta) == int: + self.save.root._v_attrs.beta = self.pp_params.disk.beta else: - print("ERROR: {} is unknown".format(name)) - self.save.close() - - def _process_rule(self, name, rule, ax_los, overwrite): - name_full = rule.group + '/' + name + '_' + ax_los - if rule.is_valid(self.save, ax_los): - if overwrite or not name_full in self.save: - data = rule.process(ax_los) - self._save_data(name_full, data, rule.description) - else: - print("Data for {} is already computed, skipping...".format(name_full)) - else: - print("ERROR: {} is not valid in this context".format(name_full)) - - - def _save_data(self, name_full, data, description): - """ - Save data in the HDF5 structure, overwrite if necessary - """ - if name_full in self.save: - node = self.save.get_node(name_full) - del node - self.save.create_array(os.path.dirname(name_full), os.path.basename(name_full), - data, description, createparents=True) + self.save.root._v_attrs.beta = int(self.save.root._v_attrs.run.split('_')[1][4:]) def _coldens(self, ax_los): datamap = self._rt.process(self._cam[ax_los], surf_qty=True) @@ -211,6 +296,12 @@ class PostProcessor: dmap_jeans_ratio = dmap_jeans * 2**(dmap_levels) return dmap_jeans_ratio + def _jeans_ratio(self, ax_los): + dmap_jeans = self.save.get_node('/maps/jeans_' + ax_los).read() + dmap_levels = self.save.get_node('/maps/levels_' + ax_los).read() + dmap_jeans_ratio = dmap_jeans * 2**(dmap_levels) + return dmap_jeans_ratio + def _toomreQ_disk(self, ax_los): """ Compute the Toomre Q parameter in a Keplerian disk @@ -256,21 +347,338 @@ class PostProcessor: return map_Q + def _radial_bins(self, ax_los): + pos_star = self.pp_params.disk.pos_star + im_extent = self.save.root.maps._v_attrs.im_extent + + # radius of the corner of the box plus a margin + rad_of_box = np.sqrt((im_extent[1] - pos_star[0])**2 + (im_extent[3] - pos_star[1])**2) + 0.1 + + bin_in = self.pp_params.disk.bin_in + bin_out = self.pp_params.disk.bin_out + nb_bin = self.pp_params.disk.nb_bin + + # radial bins + if self.pp_params.disk.binning == 'log': + lrad_in = np.log10(bin_in) + lrad_ext = np.log10(bin_out) + rad_bins = np.logspace(lrad_in, lrad_ext, num=nb_bin) + elif binning == 'lin': + rad_bins = np.linspace(bin_in, bin_out, num=nb_bin) + + # Add boundaries + rad_bins = np.concatenate(([0.], rad_bins, [rad_of_box])) + return rad_bins + + def _rr(self, ax_los): + im_extent = self.save.root.maps._v_attrs.im_extent + map_size = self.pp_params.out.map_size + pos_star = self.pp_params.disk.pos_star + + x = np.linspace(im_extent[0], im_extent[1], map_size) + y = np.linspace(im_extent[2], im_extent[3], map_size) + xx, yy = np.meshgrid(x, y) + rr = np.sqrt((xx - pos_star[0])**2 + (yy - pos_star[1])**2) + return rr + + def _bins_on_map(self, ax_los): + rad_bins = self.save.get_node('/radial/radial_bins_' + ax_los).read() + rr = self.save.get_node('/maps/rr_' + ax_los).read() + + # Find appropriate bin for each coordinate set + bins = np.zeros(rr.shape, dtype=int) + for r in rad_bins[1:]: + bins = bins + (rr >= r).astype(int) + return bins + + def _rad_avg(self, name, ax_los): + radial_bins = self.save.get_node('/radial/radial_bins_' + ax_los).read() + bins_on_map = self.save.get_node('/maps/bins_on_map_' + ax_los).read() + dmap = self.save.get_node('/maps/' + name + '_' + ax_los).read() + + # mean of all the cells in the bin + mean_bin = np.zeros(len(radial_bins) - 1) + for j in range(len(radial_bins) - 1): + mean_bin[j] = np.mean(dmap[bins_on_map == j]) + return mean_bin + + def _rad_avg_map(self, name, ax_los): + + radial_bins = self.save.get_node('/radial/radial_bins_' + ax_los).read() + bins_on_map = self.save.get_node('/maps/bins_on_map_' + ax_los).read() + rr = self.save.get_node('/maps/rr_' + ax_los).read() + mean_bin = self.save.get_node('/radial/rad_avg_' + name + '_' + ax_los).read() + + # Add value for border + mean_bin = np.concatenate(([mean_bin[0]], mean_bin)) + + rr_flat = rr.flatten() + bins_on_map_flat = bins_on_map.flatten() + + + # Compute the map azimuthally averaged + # use linear interpolation to improve accuracy + avg_flat = (radial_bins[bins_on_map_flat + 1] - rr_flat) * mean_bin[bins_on_map_flat] + avg_flat = avg_flat + (rr_flat - radial_bins[bins_on_map_flat]) * mean_bin[bins_on_map_flat + 1] + avg_flat = avg_flat / (radial_bins[bins_on_map_flat + 1] - radial_bins[bins_on_map_flat]) + avg_map = np.reshape(avg_flat, rr.shape) + + return avg_map + + def _fluct_map(self, name, ax_los): + + dmap = self.save.get_node('/maps/' + name + '_' + ax_los).read() + avg_map = self.save.get_node('/maps/avg_map_' + name + '_' + ax_los).read() + + return dmap / avg_map + + def _pdf(self, name, ax_los): + fluct_map = self.save.get_node('/maps/fluct_' + name + '_' + ax_los).read() + rr = self.save.get_node('/maps/rr_' + ax_los).read() + + mask_pdf = (rr > self.pp_params.disk.rmin_pdf) & (rr < self.pp_params.disk.rmax_pdf) + + nb_cells = np.sum(mask_pdf.flatten()) + values, edges = np.histogram(np.log10(fluct_map[mask_pdf].flatten()), + self.pp_params.pdf.nb_bin, + weights = np.ones(nb_cells) / nb_cells) + centers = 0.5 * (edges[1:] + edges[:-1]) + return np.stack([values, centers]) + + def _fit_pdf(self, name, ax_los): + pdf = self.save.get_node('/hist/pdf_' + name + '_' + ax_los) + values, centers = pdf.read() + mask_fit = ((centers > self.pp_params.pdf.xmin_fit) & + (centers < self.pp_params.pdf.xmax_fit) & + (values > 0)) + (slope, origin, correlation, _, stderr) = linregress(centers[mask_fit], np.log10(values[mask_fit])) + + pdf.attrs.slope = slope + pdf.attrs.origin = origin + pdf.attrs.correlation = correlation + pdf.attrs.stderr = stderr + pdf.attrs.var = np.var + return True + + def _sinks(self): + csv_name = self.path + '/output_' + str(self.num).zfill(5) + '/sink_' + str(self.num).zfill(5) + '.csv' + sinks = np.loadtxt(csv_name, delimiter=',') + header = ['Id', 'M', 'dmf', 'x', 'y', 'z', 'vx', 'vy', 'vz', + 'rot_period', 'lx', 'ly', 'lz', + 'acc_rate', 'acc_lum', 'age', 'int_lum', 'Teff'] + if len(sinks) == 0: + sinks = np.zeros(len(header)) + + sinks_dict = dict() + for key, a in zip(header, sinks): + sinks_dict[key] = a + + return sinks_dict + def def_rules(self): self.rules = { - 'coldens' : Rule(self._coldens, "Column density", '/maps'), - 'rho' : Rule(self._rho, "Density slice", '/maps'), - 'speed_h' : Rule(self._speed_h, "Horizontal speed slice wrt the line of sight", '/maps'), - 'speed_v' : Rule(self._speed_v, "Vertical speed slice wrt the line of sight", '/maps'), - 'T' : Rule(self._temperature, "Temperature slice", '/maps', dependencies=['rho']), - 'levels' : Rule(self._levels, "Max level within line of sight", '/maps'), - 'jeans' : Rule(self._jeans, "Jeans lenght slice", '/maps', dependencies=['rho', 'T']), - 'jeans_ratio' : Rule(self._jeans_ratio, "Jeans' lenght divided by the max resolution", + # Base rules + 'coldens' : Rule(self, self._coldens, "Column density", '/maps'), + 'rho' : Rule(self, self._rho, "Density slice", '/maps'), + 'speed_h' : Rule(self, self._speed_h, "Horizontal speed slice wrt the line of sight", '/maps'), + 'speed_v' : Rule(self, self._speed_v, "Vertical speed slice wrt the line of sight", '/maps'), + 'T' : Rule(self, self._temperature, "Temperature slice", '/maps', dependencies=['rho']), + 'levels' : Rule(self, self._levels, "Max level within line of sight", '/maps'), + 'jeans' : Rule(self, self._jeans, "Jeans lenght slice", '/maps', dependencies=['rho', 'T']), + 'jeans_ratio' : Rule(self, self._jeans_ratio, "Jeans' lenght divided by the max resolution", '/maps', dependencies=['jeans', 'levels']), - 'Q' : Rule(self._toomreQ_disk, "Toomre Q parameter for a Keplerian disk", '/maps', - dependencies=['coldens'], axes=['z'], - is_valid=lambda save, axe: self.pp_params.disk.on) + 'Q' : Rule(self, self._toomreQ_disk, "Toomre Q parameter for a Keplerian disk", '/maps', + dependencies=['coldens'], args_ok=['z'], + is_valid=lambda _: self.pp_params.disk.on), + 'sinks' : Rule(self, self._sinks, group="/datasets", args_ok=[None], + description={'Id': '', 'M':'[Msol]', 'dmf':'[Msol]', + 'x': '', 'y': '', 'z': '', 'vx': '', 'vy': '', 'vz': '', + 'rot_period':'[y]', 'lx':'|l|', 'ly':'|l|', 'lz':'|l|', + 'acc_rate':'[Msol/y]', 'acc_lum':'[Lsol]', 'age':'[y]', + 'int_lum':'[Lsol]', 'Teff':'[K]'}), + # Helpers + 'radial_bins' : Rule(self, self._radial_bins, "Radial bins", '/radial', args_ok=['z']), + 'rr' : Rule(self, self._rr, "Coordinate map", '/maps', args_ok=['z']), + 'bins_on_map' : Rule(self, self._bins_on_map, "Convert map coordinates to bins", '/maps', + dependencies=['radial_bins', 'rr'], args_ok=['z']) + } + + # Average and other + averageables = ['coldens', 'rho', 'T', 'Q'] + for name in averageables: + self.rules['rad_avg_' + name] = Rule(self, partial(self._rad_avg, name), + "Azimuthal average of {}".format(name), '/radial', + dependencies=['radial_bins', 'bins_on_map', name], + args_ok=['z']) + + self.rules['avg_map_' + name] = Rule(self, partial(self._rad_avg_map, name), + "Interpolated map of azimuthal average of {}".format(name), + '/maps', + dependencies=['radial_bins', 'bins_on_map', + 'rr', 'rad_avg_' + name], + args_ok=['z']) + self.rules['fluct_' + name] = Rule(self, partial(self._fluct_map, name), + "Fluctuation wrt to average of {}".format(name), + '/maps', + dependencies=[name, 'avg_map_' + name], + args_ok=['z']) + self.rules['pdf_' + name] = Rule(self, partial(self._pdf, name), + "Probability density function of {} fluctuations".format(name), + '/hist', + dependencies=['rr', 'fluct_' + name], + args_ok=['z']) + + self.rules['fit_pdf_' + name] = Rule(self, partial(self._fit_pdf, name), + "Fit the PDF of {} fluctuations".format(name), + '/hist', + dependencies=['pdf_' + name], + args_ok=['z']) + +class Comparator(BaseProcessor): + """ + Do comparaison between outputs and runs + """ + + def __init__(self, path, runs, nums, path_out=None, pp_params=Params()): + """ + Creates the basic structures needed for the outputs + """ + + self.pp_params = pp_params + + # Determining output directory + if (path_out is None): + path_out = path + + # Open outfile + if not pp_params.out.tag == '': + tag_name = '_' + pp_params.out.tag + else : + tag_name = '' + + self.filename = (path_out + '/comp' + tag_name + '.h5') + self.save = tables.open_file(self.filename, mode="a", title="Comparaison file") + + # Get postprocesor objets for each run + self.pp_runs = dict() + if not type(nums) == dict: + nums_tmp = nums + nums = dict() + for run in runs: + nums[run] = nums_tmp + + for run in runs: + path_run = path + '/' + run + path_out_run = path_out + '/' + run + self.pp_runs[run] = dict() + for num in nums[run]: + self.pp_runs[run][num] = PostProcessor(path_run, num, path_out=path_out_run, pp_params=pp_params) + + # save metadata + self.save.root._v_attrs.runs = runs + self.save.root._v_attrs.nums = nums + + # log info + self.log_id = "[comp {}] ".format(self.pp_params.out.tag) + + self.save.close() + self.def_rules() + + def _time_series(self, name, getter): + nums = self.save.root._v_attrs.nums + series = dict() + for run in self.save.root._v_attrs.runs: + series[run] = np.zeros(len(nums[run])) + for i, num in enumerate(nums[run]): + series[run][i] = getter(self.pp_runs[run][num]) + return series + + def _comp(self, name, getter): + runs = self.save.root._v_attrs.runs + nums = self.save.root._v_attrs.nums + prop = np.zeros(len(runs)) + for i, run in enumerate(runs): + num = nums[run][0] + prop[i] = getter(self.pp_runs[run][num]) + return prop + + def _time_avg(self, name): + runs = self.save.root._v_attrs.runs + mean = np.zeros(len(runs)) + std = np.zeros(len(runs)) + for i, run in enumerate(runs): + serie = self.save.get_node('/series/' + name + '/' + run).read() + mean[i] = np.mean(serie) + std[i] = np.std(serie) + return {"mean": mean, "std": std} + + def _get_attr(self, attr_name, pp): + h5file = tables.open_file(pp.filename, "r") + attr = h5file.root._v_attrs[attr_name] + h5file.close() + return attr + + def _get_pdf_slope(self, name, pp): + pp.process(['fit_pdf_' + name], ['z'], overwrite=self.overwrite_dep) + h5file = tables.open_file(pp.filename, "r") + pdf = h5file.get_node('/hist/pdf_' + name +'_z') + slope = pdf.attrs.slope + h5file.close() + return slope + + def _get_sinks_mass(self, pp): + pp.process(['sinks'], overwrite=self.overwrite_dep) + h5file = tables.open_file(pp.filename, "r") + sinks_mass = h5file.get_node('/datasets/sinks/M').read() + h5file.close() + return np.sum(sinks_mass) + + def def_rules(self): + averageables = ['coldens', 'rho', 'T', 'Q'] + self.rules = { + 'beta' : Rule(self, lambda arg: self._comp("beta", partial(self._get_attr, 'beta')), group='/comp', + args_ok = [None]), + 'time_pdf_slope' : Rule(self, + lambda name: self._time_series("pdf_slope_" + name, + partial(self._get_pdf_slope, name)), + group='/series', args_ok = averageables), + 'time_sinks_mass' : Rule(self, partial(self._time_series, "sinks", self._get_sinks_mass), + group='/series', args_ok=[None]), + 'time' : Rule(self, partial(self._time_series, "time", partial(self._get_attr, 'time')), + group='/series', args_ok=[None]), + 'avg_pdf_slope' : Rule(self, + lambda name: self._time_avg("time_pdf_slope_" + name), + group='/comp', dependencies=['time_pdf_slope'], + args_ok=averageables, + description={"mean": "Temporal average", "std": "Standard deviation"}) } +def get_time(path, num): + """ + Return the time of the output (code units) + + Parameters + ---------- + num output number + path_out path of the pipeline output + + Returns + ------- + time the time of the output (code units) + """ + try: + f = open(path + '/output_' + str(num).zfill(5) + '/info_' + str(num).zfill(5) + '.txt') + for line in f: + ls = line.split() + if len(ls) > 1 and ls[0] == 'time': + time = float(ls[2]) + break + # ro = pymses.RamsesOutput(path, num, order='>') + # time = ro.info['time'] # time in codeunits + f.close() + return time + except IOError as e: + print(e) + return np.nan diff --git a/params.py b/pp_params.py similarity index 56% rename from params.py rename to pp_params.py index e7ecb17..d0444b3 100644 --- a/params.py +++ b/pp_params.py @@ -1,5 +1,7 @@ # coding: utf-8 +import numpy as np + class PlotParams: """ Plot parameters @@ -9,7 +11,7 @@ class PlotParams: ntick = 6 # Number of ticks for maps set_lim = True # Set default limits vel_red = 40 # Take point each vel_red for velocities - + put_title = False class DiskParams: @@ -18,6 +20,24 @@ class DiskParams: """ on = False # Enable specific disk analysis pos_star = np.array([1., 1., 1.]) # Position of the central star + binning = "log" # Kind of binning (lin = linear, log = logarithmic) + nb_bin = 100 # Number of bins for averaged quantities + bin_in = 1e-3 # Outer radius of the inner bin + bin_out = 0.25 # Inner radius of the outer bin + rmin_pdf = 0.075 # Inner radius for PDF computation + rmax_pdf = 0.3 # Outer radius for PDF computation + + beta = False # Beta cooling. Do nothing if False. + # If true, beta will be parsed, + # otherwise the value is read therre + +class PdfParams: + """ + parameters for probability density functions + """ + nb_bin = 50 # Number of bins for the PDF + xmin_fit = 0. # Lower boundary of the fit + xmax_fit = 1.25 # Upper boundary of the fit class PymsesParams: @@ -37,12 +57,12 @@ class OutputParams: tag = "" # Tag for the image - class Params: """ Strutured parameters for the post processing """ disk = DiskParams() + pdf = PdfParams pymses = PymsesParams() out = OutputParams() plot = PlotParams()