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

Began the work on a new pipeline using the new HDF5 version

Browse Source
This commit is contained in:
Noe Brucy
2019-11-04 16:05:09 +01:00
parent 4ff3e9f605
commit ed55a0cba1
4 changed files with 1083 additions and 220 deletions
Download Patch File Download Diff File Expand all files Collapse all files
pipeline.py
+284
View File
@@ -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)
plotter.py
+234 -83
View File
@@ -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,103 +79,114 @@ 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
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, 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):
self._plot_map(name, ax_los)
rule.plot(arg)
P.tight_layout(pad=1)
P.savefig(plot_filename)
P.close()
self._log("{} plotted".format(name_full), "SUCCESS")
else:
print("Data for {} is already computed, skipping...".format(name_full))
self._log("Plot {} is already done, skipping...".format(name_full))
else:
self._log("{} is not valid in this context".format(name_full), "ERROR")
self._file_out.close()
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',
cmap=cmap,
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.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
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]
@@ -167,25 +197,146 @@ class Plotter:
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,
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:
cbar.set_label(name)
xlabel = name_x
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)
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:
y = node_y.mean.read()
yerr = node_y.std.read()
P.errorbar(x, y, yerr=y, fmt='*')
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:
postprocessor.py
+472 -64
View File
@@ -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,12 +152,12 @@ 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
"""
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
@@ -62,16 +167,19 @@ class PostProcessor:
# Open outfile
if not pp_params.out.tag == '':
tag_name = '_' + 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'))
title=os.path.basename(path)+ '_' + format(num,'05'))
# 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"])
@@ -97,10 +205,14 @@ class PostProcessor:
time = self._ro.info['time'] # time in codeunits
# 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
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
@@ -119,52 +231,25 @@ class PostProcessor:
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.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)
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)
# 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("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
params.py → pp_params.py
+22 -2
View File
@@ -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()