From ea6f9b6bdd3a1d1dd127def98aad72a2ae343c95 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Wed, 13 Nov 2019 17:33:15 +0100 Subject: [PATCH] Added extror from logs and namelist --- pipeline.py | 37 +++-- plotter.py | 342 ++++++++++++++++++++++++++++++++-------------- postprocessor.py | 347 +++++++++++++++++++++++++++++++++++++---------- pp_params.py | 86 ++++-------- 4 files changed, 563 insertions(+), 249 deletions(-) diff --git a/pipeline.py b/pipeline.py index 3f0a2ae..67dc8c4 100644 --- a/pipeline.py +++ b/pipeline.py @@ -27,6 +27,9 @@ input_args.add_argument("-p", "--project", help="specify project name (directory within the input directory)", default="disk") +input_args.add_argument("-c", "--config", + help="Path of a default config file", + default=None) input_args.add_argument("-wo", "--which_inputs", choices=['all', 'id', 'time'], @@ -77,15 +80,19 @@ output_args.add_argument("--tag", help="Add a special tag on output filemanes", default='') +output_args.add_argument("--interactive", + help="Interactive mode : do not save data", + action='store_true') + 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) + default=False) + pp_args = parser.add_argument_group('postproc', "Post Processing configuration") @@ -118,7 +125,7 @@ pp_args.add_argument("--plot", pp_args.add_argument("-plargs", "--plot_args", help="Args to give to plot rules", - default=['x', 'y', 'z'], + default=[None], nargs='*') pp_args.add_argument("-d", "--disk", @@ -174,11 +181,15 @@ runs = args.runs storage_in = args.input_path storage_out = args.output_path -pp_params = Params() +if args.config is None: + pp_params = default_params() +else: + pp_params = load_params(args.config) pp_params.out.zoom = args.zoom pp_params.out.tag = args.tag pp_params.out.map_size = args.map_size +pp_params.out.interactive = args.interactive pp_params.pymses.fft = args.fft @@ -208,7 +219,6 @@ 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 @@ -249,15 +259,11 @@ for run in runs: while not success: try: - if len(args.process) > 0 and len(args.plot) > 0: + if len(args.process) > 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) @@ -274,10 +280,15 @@ for run in runs: else: raise -if len(args.compare) > 0: - path_in = storage_in + project - path_out = storage_out + project +path_in = storage_in + project +path_out = storage_out + project +if len(args.plot) > 0: + pl = Plotter(path_in, runs, nums_success, path_out=path_out, pp_params=pp_params) + pl.process(args.plot, args.plot_args, + overwrite=args.overwrite, overwrite_dep=args.overwrite_dependencies) + +if len(args.compare) > 0: 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 091a97b..58514ff 100644 --- a/plotter.py +++ b/plotter.py @@ -15,9 +15,10 @@ import matplotlib.patches as mpatches from matplotlib.collections import PatchCollection from functools import partial from numpy.polynomial.polynomial import polyfit +from scipy.ndimage.filters import gaussian_filter1d from pp_params import * -from postprocessor import Rule, BaseProcessor +from postprocessor import * P.rcParams['image.cmap']='plasma' @@ -29,18 +30,12 @@ P.rcParams.update(tex_params) class PlotRule(Rule): - def plot(self, arg): - return self.process_fn(arg) + def plot(self, save, arg, **kwargs): + self.postproc.save = save + return self.process_fn(arg, **kwargs) 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) + return arg in self.args_ok and self.is_valid_add(arg) class Plotter(BaseProcessor): """ @@ -55,81 +50,139 @@ class Plotter(BaseProcessor): G = 1. # Gravitational constant - def __init__(self, filename=None, path_out='.', num=None, pp_params=Params()): + def __init__(self, path, runs, nums, path_out=None, pp_params=default_params(), tag=None): - self.pp_params = pp_params + if type(pp_params) == str: + self.pp_params = load_params(pp_params) + else : + self.pp_params = pp_params + + if tag is not None: + self.pp_params.out.tag = tag # Determining output directory - if path_out is None: - if filename is None: - self.path_out='.' - else: - self.path_out = os.path.dirname(filename) + if (path_out is None): + self.path_out = path else: self.path_out = path_out - # Find HDF5 file - if filename is None: - if not pp_params.out.tag == '': - tag_name = '_' + pp_params.out.tag - else : - tag_name = '' - self.filename = (self.path_out + '/postproc_' + - tag_name + format(num,'05') + '.h5') - else: - self.filename = filename + # 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_run, pp_params) + + # Get comparator object + self.comp = Comparator(path, runs, nums, path_out, pp_params) + + # save infos + self.runs = runs + self.nums = nums + self.log_id = "[plot {}] ".format(self.pp_params.out.tag) self.def_rules() - 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) + def _process_single(self, name, rule, arg, overwrite=False, just_done=[], **kwargs): + # Solve dependencies + for dep in rule.dependencies: + if dep in self.comp.rules: + if type(rule.dependencies) == dict: + self.comp.process([dep], [rule.dependencies[dep]], + self.overwrite_dep) + else: + self.comp.process([dep], [arg], self.overwrite_dep, self.overwrite_dep) + else: + for run in self.runs: + for i, num in enumerate(self.nums[run]): + if type(rule.dependencies) == dict: + dep_arg = rule.dependencies[dep] + self.pp_runs[run][num].process([dep], [dep_arg], + self.overwrite_dep) + else: + self.pp_runs[run][num].process([dep], [arg], + self.overwrite_dep, + self.overwrite_dep) - def _process_single(self, name, rule, arg, overwrite=False, overwrite_dep=False, just_done=[]): - done = self._plot_rule(name, rule, arg, overwrite) - return [] + # Process rule + done = self._process_rule(name, rule, arg, overwrite, just_done, **kwargs) + return just_done + [done] - def _plot_rule(self, name, rule, arg, overwrite): + def _process_rule(self, name, rule, arg, overwrite=False, just_done=[], **kwargs): if not arg is None: - name_full = rule.group + '/' + name + '_' + str(arg) + name_full = name + '_' + str(arg) else: - name_full = rule.group + '/' + name + name_full = 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") + if rule.kind == 'classic': + for run in self.runs: + for i, num in enumerate(self.nums[run]): + plot_filename = self._find_filename(name_full, run, num) + save = tables.open_file(self.pp_runs[run][num].filename) + try: + self._plot_rule(rule, save, arg, plot_filename, overwrite, **kwargs) + finally: + save.close() else: - self._log("Plot {} is already done, skipping...".format(name_full)) + if rule.kind == 'series' and len(self.runs) == 1: + run = self.runs[0] + plot_filename = self._find_filename(name_full, run) + else: + plot_filename = self._find_filename(name_full) + save = tables.open_file(self.comp.filename, 'r') + try : + self._plot_rule(rule, save, arg, plot_filename, overwrite, **kwargs) + finally: + save.close() else: self._log("{} is not valid in this context".format(name_full), "ERROR") + def _plot_rule(self, rule, save, arg, plot_filename, overwrite, **kwargs): + if overwrite or not os.path.exists(plot_filename): + if not self.pp_params.out.interactive: + P.figure() + rule.plot(save, arg, **kwargs) + P.tight_layout(pad=1) + if not self.pp_params.out.interactive: + P.savefig(plot_filename) + P.close() + self._log("{} plotted".format(plot_filename), "SUCCESS") + else: + self._log("{} plotted".format(os.path.basename(plot_filename)), "SUCCESS") + else: + self._log("Plot {} is already done, skipping...".format(plot_filename)) - def _find_filename(self, name_full): + + def _find_filename(self, name_full, run=None, num=None): 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 + + if not run is None and not num is None: + return (self.path_out + '/' + run + '/' + name_full + tag_name + '_' + format(num,'05') + self.pp_params.plot.out_ext) + elif not run is None: + return (self.path_out + '/' + run + '/' + + name_full + tag_name + 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, 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.save.root.maps._v_attrs.im_extent @@ -201,8 +254,8 @@ class Plotter(BaseProcessor): P.quiverkey(Q, 0.6, 0.98, max_v, r'$'+str(max_v)[0:4]+'$ (code)', labelpos='E', coordinates='figure') - def _plot_radial(self, name, ax_los='z', label=None, xlog=False, ylog=False): - P.figure() + def _plot_radial(self, name, label=None, xlog=False, ylog=False): + 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() @@ -220,7 +273,7 @@ class Plotter(BaseProcessor): 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] @@ -236,106 +289,185 @@ class Plotter(BaseProcessor): P.ylim([None, 1.]) - def _plot(self, name_x, name_y, xlabel=None, ylabel=None, linearfit=False): - P.figure() + def _plot(self, name_x, name_y, xlabel=None, ylabel=None, + xunit=None, yunit=None, + linearfit=False, smooth=0, + nml_key=None, runs=None, **kwargs): + 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 + if 'label' in node_x._v_attrs: + xlabel = node_x._v_attrs.label else: - xlabel = name_x + xlabel = os.path.basename(name_x) if ylabel is None: - if 'label' in node_y: - ylabel = name_y._vattrs.label + if 'label' in node_y._v_attrs: + ylabel = node_y._v_attrs.label else: - ylabel = name_y + ylabel = os.path.basename(name_y) - x = node_x.read() - if node_y._v_attrs.CLASS == 'ARRAY': - y = node_y.read() - P.plot(x, y, fmt='*') + if 'unit' in node_x._v_attrs: + xunit_old = node_x._v_attrs.unit else: - y = node_y.mean.read() - yerr = node_y.std.read() - P.errorbar(x, y, yerr=y, fmt='*') + xunit_old = cst.none + + if 'unit' in node_y._v_attrs: + yunit_old = node_y._v_attrs.unit + else: + yunit_old = cst.none + + if xunit is None: + xunit = xunit_old + if yunit is None: + yunit = yunit_old + + xlabel = xlabel + unit_str(xunit) + ylabel = ylabel + unit_str(yunit) P.xlabel(xlabel) P.ylabel(ylabel) + P.grid() + + if node_y._v_attrs.CLASS == 'ARRAY': + x = node_x.read() * xunit_old.express(xunit) + y = node_y.read() * yunit_old.express(yunit) + if smooth > 0: + y = gaussian_filter1d(y, sigma=smooth) + yerr = None + P.plot(x, y, fmt='*', **kwargs) + elif 'mean' in node_y: + x = node_x.read() * xunit_old.express(xunit) + y = node_y.mean.read() * yunit_old.express(yunit) + if smooth > 0: + y = gaussian_filter1d(y, sigma=smooth) + yerr = node_y.std.read() * yunit_old.express(yunit) + P.errorbar(x, y, yerr=yerr, fmt='*', **kwargs) + else: + yerr = None + if runs is None: + runs = self.runs + for run in runs: + x_run, y_run = node_x[run], node_y[run] + x = x_run.read() * xunit_old.express(xunit) + y = y_run.read() * yunit_old.express(yunit) + if smooth > 0: + y = gaussian_filter1d(y, sigma=smooth) + if nml_key is None: + label_run = r"{}".format(self.save.root._v_attrs.attrs[y_run.name].label) + else: + prop_name = os.path.basename(nml_key) + prop_value = self.comp.get_nml(nml_key, run) + label_run = r"{} = {}".format(prop_name, prop_value) + P.plot(x, y, label=label_run, **kwargs) + P.legend() if linearfit: - if node_y._v_attrs.CLASS == 'ARRAY': - (a, b, rho, _, stderr) = linregress(node_x.read(), node_y.read()) + _overlay_linearfit(x, y, yerr) + + def _overlay_linearfit(x, y, yerr=None): + if yerr is None: + (a, b, rho, _, stderr) = linregress(x, y) else: - c = polyfit(node_x.read(), node_y.mean.read(), 1, - w = [1.0 / ty for ty in node_y.std.read()], full=False) + c = polyfit(x, y, 1, + w = [1.0 / ty for ty in yerr], 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']), + self._plot_map('coldens', ax_los, + label=r'$\Sigma$ (code)', vmin=0.01, vmax=100), + "Column density", dependencies=['coldens']), 'rho' : PlotRule(self, lambda ax_los: self._plot_map('rho', ax_los, label=r'$\rho$ (code)'), - "Density slice", ['/maps/rho']), + "Density slice", dependencies=['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']), + self._plot_map('coldens', ax_los, + label=r'$\Sigma$ (code)', + vmin=0.01, vmax=100, + overlays=[self._overlay_levels]), + "Column density", + dependencies=['coldens', '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']), + self._plot_map('rho', ax_los, label=r'$\rho$ (code)', + overlays=[self._overlay_speed]), + "Density slice", + dependencies= ['rho', 'speed_h', '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']), + dependencies=['jeans_ratio', 'levels']), 'Q' : PlotRule(self, lambda ax_los: - self._plot_map('rho', ax_los, label=r'$Q$', vmin=0.01, vmax=100, cmap='RdBu_r'), + 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']) + dependencies=['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), + 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], + dependencies=['radial_bins', '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['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=['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], + dependencies=['pdf_' + name], args_ok=['z']) self.rules.update({ - 'kappa_beta' : PlotRule(self, partial(self._plot, '/comp/beta', '/comp/avg_pdf_slope_coldens', + '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']) + args_ok=[None], kind='comp', + dependencies=['beta', 'avg_pdf_slope_coldens']), + 'sink_mass' : PlotRule(self, partial(self._plot, + '/series/sinks_from_log/time', + '/series/sinks_from_log/mass_sink', + ylabel="Mass of sinks", + xunit=cst.Myr, + yunit=cst.Msun), + args_ok=[None], kind='series', + dependencies=['sinks_from_log']), + 'assfr' : PlotRule(self, partial(self._plot, + '/series/sfr_from_log/time', + '/series/sfr_from_log/sfr', + ylabel="Averaged surfacic SFR", + xunit=cst.Myr, + yunit=cst.ssfr), + args_ok=[None], kind='series', + dependencies=['sfr_from_log']), + 'issfr' : PlotRule(self, partial(self._plot, + '/series/sinks_from_log/time', + '/series/sinks_from_log/issfr', + ylabel="Surfacic SFR", + xunit=cst.Myr, + yunit=cst.ssfr), + args_ok=[None], kind='series', + dependencies=['issfr']) }) diff --git a/postprocessor.py b/postprocessor.py index a5b12dd..3c9df57 100644 --- a/postprocessor.py +++ b/postprocessor.py @@ -2,25 +2,48 @@ import sys import os +import glob as glob + import tables import pymses import numpy as np from numpy.polynomial.polynomial import polyfit +from scipy.stats import linregress + from pymses.sources.ramses import output 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 +import pymses.utils.constants as cst + +import f90nml +import bunch + from functools import partial from abc import ABCMeta, abstractmethod from pp_params import * +def unit_str(unit): + if unit == cst.none: + return '' + elif len(unit.name) > 0 : + return ' ({})'.format(unit.name) + else: + return ' ' + str(unit) + +cst.ssfr = cst.create_unit('Msun/yr/pc^2', + base_unit=cst.Msun/cst.year/cst.pc**2, + descr="Surfacic SFR", + latex='M_{\astrosun}.yr^{-1}.pc^{-2}') + + class Rule: - def __init__(self, postproc, process, description='', group='', dependencies=[], args_ok=['x', 'y', 'z'], - is_valid=lambda arg:True): + def __init__(self, postproc, process, description='', + group='', dependencies=[], args_ok=['x', 'y', 'z'], + is_valid=lambda arg:True, kind='classic', unit=cst.none): self.postproc = postproc self.process_fn = process self.dependencies = dependencies @@ -28,12 +51,14 @@ class Rule: self.group = group self.args_ok = args_ok self.description = description + self.unit=unit + self.kind = kind - def process(self, arg): + def process(self, arg, **kwargs): if not arg is None: - return self.process_fn(arg) + return self.process_fn(arg, **kwargs) else: - return self.process_fn() + return self.process_fn(**kwargs) def is_valid(self, arg): save = self.postproc.save @@ -55,7 +80,9 @@ class BaseProcessor: log_id = "" rules = dict() + data = dict() filename = "" + save = None @abstractmethod def __init__(self): @@ -67,39 +94,45 @@ class BaseProcessor: else: print(self.log_id + string) - def process(self, to_process_list, args=[None], overwrite=False, overwrite_dep=None): + def open(self): + self.save = tables.open_file(self.filename, mode="a") + + def close(self): + self.save.close() + + def process(self, to_process_list, args=[None], overwrite=False, overwrite_dep=False, **kwargs): """ 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) + just_done = self._process_single(name, rule, arg, overwrite, just_done, **kwargs) 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=[]): + def _process_single(self, name, rule, arg, overwrite=False, just_done=[], **kwargs): # 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) + 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) + self.open() + try: + done = self._process_rule(name, rule, arg, overwrite, just_done, **kwargs) + finally: + self.close() return just_done + [done] - def _process_rule(self, name, rule, arg, overwrite=False, just_done=[]): + def _process_rule(self, name, rule, arg, overwrite=False, just_done=[], **kwargs): if not arg is None: name_full = rule.group + '/' + name + '_' + str(arg) else: @@ -107,10 +140,11 @@ class BaseProcessor: if rule.is_valid(arg): if not name_full in just_done: - if overwrite or not name_full in self.save: + if (overwrite or + (not name_full in self.save and not name_full in self.data)): self._log("Processing {}".format(name_full)) - data = rule.process(arg) - self._save_data(name_full, data, rule.description) + data = rule.process(arg, **kwargs) + self._save_data(name_full, data, rule.description, rule.unit) self._log("Data for {} computed".format(name_full), "SUCCESS") return name_full else: @@ -119,24 +153,43 @@ class BaseProcessor: self._log("{} is not valid in this context".format(name_full), "ERROR") - def _save_data(self, name_full, data, description): + def _save_data(self, name_full, data, description, unit): """ 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: + if type(data) == dict: + if type(description) == str: + self.save.create_group(os.path.dirname(name_full), + os.path.basename(name_full), + description, createparents=True) + else: + self.save.create_group(os.path.dirname(name_full), + os.path.basename(name_full), + '', createparents=True) + + if not type(unit) == dict: + self.save.get_node(name_full)._v_attrs.unit = unit + for key in data: if type(description) == dict: - self.save.create_array(name_full, key, - data[key], description[key], createparents=True) + if type(unit) == dict: + self._save_data(name_full + '/' + key, data[key], description[key], unit[key]) + else: + self._save_data(name_full + '/' + key, data[key], description[key], unit) else: - self.save.create_array(name_full, key, - data[key], description, createparents=True) + if type(unit) == dict: + self._save_data(name_full + '/' + key, data[key], '', unit[key]) + else: + self._save_data(name_full + '/' + key, data[key], '', unit) + else: + self.save.create_array(os.path.dirname(name_full), + os.path.basename(name_full), + data, description, createparents=True) + self.save.get_node(name_full).attrs.unit = unit + @abstractmethod def def_rules(self): @@ -154,14 +207,18 @@ class PostProcessor(BaseProcessor): G = 1. # Gravitational constant - def __init__(self, path=None, num=None, path_out=None, pp_params=Params()): + def __init__(self, path=None, num=None, path_out=None, pp_params=default_params()): """ Creates the basic structures needed for the outputs """ 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 + + if type(pp_params) == str: + self.pp_params = load_params(pp_params) + else : + self.pp_params = pp_params # Determining output directory if (path_out is None): @@ -175,6 +232,9 @@ class PostProcessor(BaseProcessor): self.filename = (path_out + '/postproc_' + tag_name + format(num,'05') + '.h5') + + if not os.path.exists(path_out): + os.makedirs(path_out) self.save = tables.open_file(self.filename, mode="a", title=os.path.basename(path)+ '_' + format(num,'05')) @@ -184,6 +244,13 @@ class PostProcessor(BaseProcessor): self.num = num self._ro = pymses.RamsesOutput(path, num, order=pp_params.pymses.order) self._amr = self._ro.amr_source(["rho","vel","P"]) + self.info = self._ro.info.copy() + + # Delete specific info + self.info.pop('dom_decomp_Hilbert_keys') + self.info.pop('nstep_coarse') + self.info.pop('dom_decomp') + self.info.pop('time') # Density operator self._rho_op = ScalarOperator(lambda dset: dset["rho"], self._ro.info["unit_density"]) @@ -196,7 +263,7 @@ class PostProcessor(BaseProcessor): # Set the extend of the image self._radius = 0.5 / pp_params.out.zoom - self._lbox = self._ro.info['boxlen'] + self._lbox = self.info['boxlen'] center = pp_params.out.center im_extent = [(- self._radius + center[0]) * self._lbox, ( self._radius + center[0]) * self._lbox, @@ -219,7 +286,6 @@ class PostProcessor(BaseProcessor): 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 @@ -246,12 +312,15 @@ class PostProcessor(BaseProcessor): Add additional metadata to the file """ - # 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: - self.save.root._v_attrs.beta = int(self.save.root._v_attrs.run.split('_')[1][4:]) + # Label of the run in the label.txt file + label_filename = self.path + '/' + self.pp_params.input.label_filename + if os.path.exists(label_filename): + label_file = open(label_filename, 'r') + self.label = label_file.readline() + label_file.close() + else: + self.label = self.run + self.save.root._v_attrs.label = self.label def _coldens(self, ax_los): datamap = self._rt.process(self._cam[ax_los], surf_qty=True) @@ -482,27 +551,39 @@ class PostProcessor(BaseProcessor): # 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", + '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, self._toomreQ_disk, "Toomre Q parameter for a Keplerian disk", '/maps', + '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]', + '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']), + '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', + 'bins_on_map' : Rule(self, self._bins_on_map, + "Convert map coordinates to bins", '/maps', dependencies=['radial_bins', 'rr'], args_ok=['z']) } @@ -542,12 +623,15 @@ class Comparator(BaseProcessor): Do comparaison between outputs and runs """ - def __init__(self, path, runs, nums, path_out=None, pp_params=Params()): + def __init__(self, path, runs, nums, path_out=None, pp_params=default_params()): """ Creates the basic structures needed for the outputs """ - self.pp_params = pp_params + if type(pp_params) == str: + self.pp_params = load_params(pp_params) + else : + self.pp_params = pp_params # Determining output directory if (path_out is None): @@ -570,16 +654,31 @@ class Comparator(BaseProcessor): for run in runs: nums[run] = nums_tmp + attrs = dict() + self.namelist = dict() + 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) + self.pp_runs[run][num] = PostProcessor(path_run, num, path_out=path_out_run, + pp_params=pp_params) + + h5file = tables.open_file(self.pp_runs[run][nums[run][0]].filename, 'r') + attrs[run] = h5file.root._v_attrs + h5file.close() + path_nml = path_run + '/' + self.pp_params.input.nml_filename + self.namelist[run] = bunch.bunchify(f90nml.read(path_nml)) + + + self.info = self.pp_runs[runs[0]][nums[runs[0]][0]].info # save metadata self.save.root._v_attrs.runs = runs self.save.root._v_attrs.nums = nums + self.save.root._v_attrs.attrs = attrs + self.save.root._v_attrs.namelist = bunch.unbunchify(self.namelist) # log info self.log_id = "[comp {}] ".format(self.pp_params.out.tag) @@ -593,7 +692,7 @@ class Comparator(BaseProcessor): 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]) + series[run][i] = getter(run, num) return series def _comp(self, name, getter): @@ -602,7 +701,7 @@ class Comparator(BaseProcessor): prop = np.zeros(len(runs)) for i, run in enumerate(runs): num = nums[run][0] - prop[i] = getter(self.pp_runs[run][num]) + prop[i] = getter(run, num) return prop def _time_avg(self, name): @@ -615,13 +714,21 @@ class Comparator(BaseProcessor): std[i] = np.std(serie) return {"mean": mean, "std": std} - def _get_attr(self, attr_name, pp): + def get_attr(self, attr_name, run, num): + pp = self.pp_runs[run][num] 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): + def get_nml(self, nml_key, run): + res = self.namelist[run] + for key in nml_key.split('/'): + res = res[key] + return res + + def get_pdf_slope(self, name, run, num): + pp = self.pp_runs[run][num] 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') @@ -629,34 +736,136 @@ class Comparator(BaseProcessor): h5file.close() return slope - def _get_sinks_mass(self, pp): + def get_sinks_mass(self, run, num): + pp = self.pp_runs[run][num] 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 _extract_sinks_from_log(self, series, log_filename, run): + cmd_grep = "grep 'Number of sink' {} -A 2".format(log_filename) + content = os.popen(cmd_grep).readlines() + for i in range(0, len(content), 4): + series['nb_sink'][run].append(np.int(content[i].split('=')[1])) + series['mass_sink'][run].append(np.float(content[i + 1].split('=')[1])) + series['time'][run].append(np.float(content[i + 2].split('=')[1])) + return series + + def _extract_sfr_from_log(self, series, log_filename, run): + cmd_grep = "grep '\[SFR' {} ".format(log_filename) + content = os.popen(cmd_grep).readlines() + for i in range(0, len(content)): + time = np.float(content[i].split(']')[0].split('=')[1].split()[0]) + sfr = np.float(content[i].split(']')[1].split('=')[1].split()[0]) + series['time'][run].append(time) + series['sfr'][run].append(sfr) + return series + + def _from_log(self, keys, extractor): + nums = self.save.root._v_attrs.nums + + # Initialize series + series = dict() + for key in keys: + series[key] = dict() + + for run in self.save.root._v_attrs.runs: + # Initialize list + for key in keys: + series[key][run] = [] + + # get one preprocessor + num = nums[run][0] + pp = self.pp_runs[run][num] + + # Get list of run files + log_files = (pp.path + '/' + + self.pp_params.input.log_prefix + '*') + + # Parse files + for log_filename in glob.glob(log_files): + series = extractor(series, log_filename, run) + + # Numpify the lists + for key in series: + series[key][run] = np.array(series[key][run]) + + # Sort the arrays + ind_sort = series['time'][run].argsort() + for key in series: + series[key][run] = series[key][run][ind_sort] + return series + + + def _ssfr_from_mass_sink(self): + time_unit = self.save.get_node('/series/sinks_from_log/time')._v_attrs.unit + mass_unit = self.save.get_node('/series/sinks_from_log/mass_sink')._v_attrs.unit + # Surface of the box in pc^2 + surface = (self.info['unit_length'].express(cst.pc))**2 + # WARNING : We do not multiply by boxlen since already done in 'unit_length' + ssfr = dict() + for run in self.save.root._v_attrs.runs: + time = self.save.get_node('/series/sinks_from_log/time/' + run).read() + time = time * time_unit.express(cst.year) + mass_sink = self.save.get_node('/series/sinks_from_log/mass_sink/' + run).read() + mass_sink = mass_sink * mass_unit.express(cst.Msun) + sfr = (mass_sink[1:] - mass_sink[:-1]) / (time[1:] - time[:-1]) + ssfr[run] = np.zeros(len(mass_sink)) + ssfr[run][1:] = sfr / surface + return ssfr + + 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]), + # Read from log + 'sinks_from_log' : Rule(self, + partial(self._from_log, + ['time', 'mass_sink', 'nb_sink'], + self._extract_sinks_from_log), + group='/series', args_ok=[None], + unit={'time' : self.info['unit_time'], + 'mass_sink' : cst.Msun, + 'nb_sink' : cst.none}), + 'issfr' : Rule(self, + self._ssfr_from_mass_sink, + group='/series/sinks_from_log', args_ok=[None], + unit=cst.ssfr, + description="Instantaneous surfacic star formation rate", + dependencies=['sinks_from_log']), + 'sfr_from_log' : Rule(self, + partial(self._from_log, + ['time', 'sfr'], + self._extract_sfr_from_log), + group='/series', args_ok=[None], + unit={'time' : cst.year, + 'sfr' : cst.ssfr}, + description={'time': 'Time', + 'sfr' : 'Averaged surfacic star formation rate'}), + # Read from outputs + 'time' : Rule(self, partial(self._time_series, "time", + partial(self.get_attr, 'time')), + group='/series', args_ok=[None]), 'time_pdf_slope' : Rule(self, - lambda name: self._time_series("pdf_slope_" + name, - partial(self._get_pdf_slope, name)), + 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]), + # Averages and run properties 'avg_pdf_slope' : Rule(self, - lambda name: self._time_avg("time_pdf_slope_" + name), + 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"}) + description={"mean": "Temporal average", + "std": "Standard deviation"}), } + def get_time(path, num): """ Return the time of the output (code units) diff --git a/pp_params.py b/pp_params.py index d0444b3..bb3500e 100644 --- a/pp_params.py +++ b/pp_params.py @@ -1,69 +1,31 @@ # coding: utf-8 import numpy as np +import re +import bunch +import yaml +import os -class PlotParams: - """ - Plot parameters - """ - out_ext = '.jpeg' # extension for plots - put_title = False # Add a title to plot - 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 +_dir_path = os.path.dirname(os.path.realpath(__file__)) + +# Add support for '1e3' kind of float +_loader = yaml.SafeLoader +_loader.add_implicit_resolver( + u'tag:yaml.org,2002:float', + re.compile(u'''^(?: + [-+]?(?:[0-9][0-9_]*)\\.[0-9_]*(?:[eE][-+]?[0-9]+)? + |[-+]?(?:[0-9][0-9_]*)(?:[eE][-+]?[0-9]+) + |\\.[0-9_]+(?:[eE][-+][0-9]+)? + |[-+]?[0-9][0-9_]*(?::[0-5]?[0-9])+\\.[0-9_]* + |[-+]?\\.(?:inf|Inf|INF) + |\\.(?:nan|NaN|NAN))$''', re.X), + list(u'-+0123456789.')) -class DiskParams: - """ - Disk speficic parameters - """ - 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: - """ - Parameters for Pymses reader - """ - order = '<' # In which order the output are read - fft = False # Quick and dirty rendering using FFT - -class OutputParams: - """ - Parameters for post processing - """ - center = [0.5, 0.5, 0.5] # Center of the image - zoom = 1. # Zoom of the image - map_size = 512 # Size of the computed maps in pixel - - tag = "" # Tag for the image - -class Params: - """ - Strutured parameters for the post processing - """ - disk = DiskParams() - pdf = PdfParams - pymses = PymsesParams() - out = OutputParams() - plot = PlotParams() +def load_params(filename): + with open(filename) as f: + para_disk = yaml.load(f, Loader=_loader) + return bunch.bunchify(para_disk) +def default_params(): + return load_params(_dir_path + '/pp_params.yml')