diff --git a/aggregator.py b/aggregator.py index 3ae8db4..8baec73 100644 --- a/aggregator.py +++ b/aggregator.py @@ -2,7 +2,7 @@ import numpy as np from functools import partial -import postprocessor +import snapshotprocessor try: from mpi4py.futures import MPIPoolExecutor @@ -14,44 +14,44 @@ except ModuleNotFoundError: mpi = False -def _map_aux(fun, path, path_out, pp_params, run_num, **kwargs): +def _map_aux(fun, path, path_out, params, run_num, **kwargs): try: - pp = postprocessor.PostProcessor( - path + "/" + run_num[0], run_num[1], path_out + "/" + run_num[0], pp_params + snap = snapshotprocessor.SnapshotProcessor( + path + "/" + run_num[0], run_num[1], path_out + "/" + run_num[0], params ) except Exception as e: print(e) raise - return fun(pp, **kwargs) + return fun(snap, **kwargs) -def _map_rule(pp, rule, arg, overwrite, overwrite_dep): - return pp.process(rule, arg, overwrite, overwrite_dep) +def _map_rule(snap, rule, arg, overwrite, overwrite_dep): + return snap.process(rule, arg, overwrite, overwrite_dep) class Aggregator: - def get_pp_list(self, select=None): + def get_snap_list(self, select=None): if select is not None: runs, nums = self.selector.select(**select) else: runs = self.runs nums = self.nums - return [self.pp[run][num] for run in runs for num in nums[run]] + return [self.snaps[run][num] for run in runs for num in nums[run]] def map(self, func, select=None, num_process=None, **kwargs): - pp_list = self.get_pp_list(select) + snaps = self.get_snap_list(select) if num_process is None: - num_process = self.pp_params.process.num_process + num_process = self.params.process.num_process if num_process == 1: - result = [func(pp, **kwargs) for pp in pp_list] + result = [func(snap, **kwargs) for snap in snaps] else: - run_num = [(pp.run, pp.num) for pp in pp_list] + run_num = [(snap.run, snap.num) for snap in snaps] map_fn = partial( - _map_aux, func, self.path, self.path_out, self.pp_params, **kwargs + _map_aux, func, self.path, self.path_out, self.params, **kwargs ) if mpi: executor = MPIPoolExecutor(max_workers=num_process) diff --git a/baseprocessor.py b/baseprocessor.py index 9f8b8b1..25522c5 100644 --- a/baseprocessor.py +++ b/baseprocessor.py @@ -11,7 +11,7 @@ import numpy as np import tables from tables import HDF5ExtError -from pp_params import default_params, load_params +from params import default_params, load_params from units import U @@ -52,16 +52,16 @@ class BaseProcessor: rules = {} solve_self_dep = True - def __init__(self, path, path_out=None, pp_params=None, tag=None): - if pp_params is None: - self.pp_params = default_params() - elif type(pp_params) == str: - self.pp_params = load_params(pp_params) + def __init__(self, path, path_out=None, params=None, tag=None): + if params is None: + self.params = default_params() + elif type(params) == str: + self.params = load_params(params) else: - self.pp_params = copy.deepcopy(pp_params) + self.params = copy.deepcopy(params) if tag is not None: - self.pp_params.out.tag = tag + self.params.out.tag = tag # Determining output directory if path_out is None: @@ -70,7 +70,7 @@ class BaseProcessor: self.path_out = path_out def _log(self, string, status=""): - if self.pp_params.process.verbose: + if self.params.process.verbose: if len(status) > 0: print(status + ": " + self.log_id + string) else: diff --git a/fragdisk_galactica.py b/fragdisk_galactica.py index 946c038..1848e9f 100644 --- a/fragdisk_galactica.py +++ b/fragdisk_galactica.py @@ -14,7 +14,7 @@ from astrophysix.simdm.datafiles import PlotInfo, PlotType from astrophysix.utils.file import FileType from plotter import Plotter -from pp_params import default_params +from params import default_params from ramses_astrophysix import ramses @@ -32,18 +32,18 @@ select = { } -pp_params = default_params() -pp_params.input.nml_filename = "disk.nml" -pp_params.pymses.map_size = 2048 -pp_params.pymses.zoom = 4 -pp_params.pymses.filter = False -pp_params.pymses.variables = ["rho", "vel", "P", "g"] -pp_params.pymses.multiprocessing = True -pp_params.process.verbose = True -pp_params.disk.enable = True -pp_params.disk.nb_bin = 100 -pp_params.pdf.nb_bin = 100 -pp_params.process.num_process = 10 +params = default_params() +params.input.nml_filename = "disk.nml" +params.pymses.map_size = 2048 +params.pymses.zoom = 4 +params.pymses.filter = False +params.pymses.variables = ["rho", "vel", "P", "g"] +params.pymses.multiprocessing = True +params.process.verbose = True +params.disk.enable = True +params.disk.nb_bin = 100 +params.pdf.nb_bin = 100 +params.process.num_process = 10 in_dir = "/drf/projets/alfven-data/nbrucy/simus/fragdisk" out_dir = "/dsm/anais/storageA/nbrucy/visus/fragdisk/mnras" in_dir_conv = "/drf/projets/alfven-data/nbrucy/simus/conv_disk" @@ -166,7 +166,7 @@ pl_units = Plotter( filter_name="104_beta4_jr13", in_nums="last", path_out=out_dir, - pp_params=pp_params, + params=params, ) info = pl_units.comp.info rd = U.Unit.create_unit( @@ -223,8 +223,8 @@ for group in groups: if group == "jr13_tic": # JR13_TIC - pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}" - pp_params.astrophysix.descr_fmt = """ + params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}" + params.astrophysix.descr_fmt = """

Group {tag:.8}, $\\beta$ = {nml[cloud_params/beta_cool]}.

""" @@ -236,14 +236,14 @@ for group in groups: sort_run_by=nml_key, path_out=out_dir, tag="jr13_tic_mnras", - pp_params=pp_params, + params=params, unit_time=orp, ) if group == "jr12_tic": # JR12_TIC - pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}" - pp_params.astrophysix.descr_fmt = """ + params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.8}" + params.astrophysix.descr_fmt = """

Group {tag:.8}, $\\beta$ = {nml[cloud_params/beta_cool]}.

""" runs_12 = "0[0-9][0-9]_beta*_jr12" @@ -255,14 +255,14 @@ for group in groups: filter_nml=("cloud_params", "!=", 7), path_out=out_dir, tag="jr12_tic_mnras", - pp_params=pp_params, + params=params, unit_time=orp, ) if group == "jr12": - pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}" - pp_params.astrophysix.descr_fmt = """ + params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}" + params.astrophysix.descr_fmt = """

Group {tag:.4}, $\\beta$ = {nml[cloud_params/beta_cool]}

""" @@ -275,14 +275,14 @@ for group in groups: sort_run_by=nml_key, path_out=out_dir, tag="jr12_mnras", - pp_params=pp_params, + params=params, unit_time=orp, ) if group == "jr11": - pp_params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}" - pp_params.astrophysix.descr_fmt = """ + params.astrophysix.simu_fmt = "beta{nml[cloud_params/beta_cool]:g}_{tag:.4}" + params.astrophysix.descr_fmt = """

Group {tag:.4}, $\\beta$ = {nml[cloud_params/beta_cool]}

""" @@ -295,7 +295,7 @@ for group in groups: filter_nml=("cloud_params/beta_cool", ">", 3), path_out=out_dir_conv, tag="jr11_mnras", - pp_params=pp_params, + params=params, unit_time=orp, ) @@ -305,8 +305,8 @@ for group in groups: # -------------------------------------------------------------------------------------------------------------------- # for pl in pls: - pl.pp_params.process.verbose = True - pl.comp.pp_params.process.verbose = True + pl.params.process.verbose = True + pl.comp.params.process.verbose = True for run in pl.runs: simu = pl.simulations[run] proj.simulations.add(simu) diff --git a/galaxy.py b/galaxy.py index a40ad4d..6135ac3 100644 --- a/galaxy.py +++ b/galaxy.py @@ -98,8 +98,8 @@ def get_coldens(pp): im_extent = np.array(pp.get_attribute("/maps", "im_extent")) im_extent *= pp.info["unit_length"].express(U.kpc) - map_size = pp.pp_params.pymses.map_size - center = np.array(pp.pp_params.disk.center) + map_size = pp.params.pymses.map_size + center = np.array(pp.params.disk.center) center *= pp.info["unit_length"].express(U.kpc) # Physical size of cells diff --git a/gui.py b/gui.py index 123561f..6561773 100644 --- a/gui.py +++ b/gui.py @@ -1,4 +1,6 @@ import matplotlib as mpl +import numpy as np +from functools import partial import matplotlib.patches as patches import matplotlib.pyplot as plt from matplotlib.lines import Line2D @@ -8,14 +10,14 @@ from matplotlib.widgets import ( CheckButtons, LassoSelector, PolygonSelector, - RadioButtons, Slider, SpanSelector, ) from scipy.stats import linregress from skimage.draw import line -from postprocessor import * +from snapshotprocessor import SnapshotProcessor +from params import default_params class DraggablePoint: @@ -184,42 +186,6 @@ class DraggableLine: self.list_points[1].line.remove() -class FitterTool: - """ - Flexible fitter tool - """ - - def __init__(self, ax, bounds=None): - self.ax = ax - self.bounds = bounds - self.bounds_selector = SpanSelector( - self.ax, self.onselect, "horizontal", useblit=True, span_stays=True - ) - self.fitline = DraggableLine(self, self.ax_gamma, [0, 0], [1, 0.3], 0.05) - - if update_map and self.gamma_3d_button.get_status()[1]: - lr = linregress(rho[rho > -4], cs[rho > -4]) - (a, b, r, _, _) = lr - self.line_gamma.clear() - del self.line_gamma - rhomin, rhomax = np.min(rho), np.max(rho) - self.line_gamma = DraggableLine( - self, - self.ax_gamma, - [rhomin, a * rhomin + b], - [rhomax, a * rhomax + b], - 0.05, - ) - print("Gamma linregress : {}".format(lr)) - else: - lps = self.line_gamma.list_points - xa, ya, xb, yb = lps[0].x, lps[0].y, lps[1].x, lps[1].y - a = (yb - ya) / (xb - xa) - - def onselect(self, vmin, vmax): - self.bounds = (vmin, vmax) - - class InteractiveGUI: """ This is a matplotlib interactive session to restrain analysis to a specific area @@ -245,7 +211,7 @@ class InteractiveGUI: self.fcs = np.copy(self.fcs_map) self.fmap[np.logical_not(self.mask)] = np.nan - ## Map + # Map plt.sca(self.ax_fluct) if first: self.im = plt.imshow( @@ -262,7 +228,7 @@ class InteractiveGUI: else: self.im.set_data(self.fmap) - ## Gamma + # Gamma plt.sca(self.ax_gamma) if self.gamma_3d_button.get_status()[0]: @@ -317,7 +283,7 @@ class InteractiveGUI: self.ax_gamma.set_title("$\Gamma$ = {:.3g}".format(self.get_gamma(a))) - ## PDF + # PDF if first or update_map: plt.sca(self.ax_pdf) nb_cells = np.sum(self.mask_flat) @@ -325,7 +291,7 @@ class InteractiveGUI: self.std_nb_cells = nb_cells values, self.edges = np.histogram( np.log10(self.fmap.flatten()[self.mask_flat]), - self.pp.pp_params.pdf.nb_bin, + self.pp.params.pdf.nb_bin, weights=np.ones(nb_cells) / self.std_nb_cells, ) edges = self.edges @@ -340,8 +306,8 @@ class InteractiveGUI: centers = 0.5 * (edges[1:] + edges[:-1]) mask_fit = ( - (centers > self.pp.pp_params.pdf.xmin_fit) - & (centers < self.pp.pp_params.pdf.xmax_fit) + (centers > self.pp.params.pdf.xmin_fit) + & (centers < self.pp.params.pdf.xmax_fit) & (values > 0) ) (a, b, r, _, _) = linregress(centers[mask_fit], np.log10(values[mask_fit])) @@ -364,21 +330,21 @@ class InteractiveGUI: self.fit.set_label(r"a = {:.3g}, $R^2$ = {:.3g}".format(a, r ** 2)) plt.legend() - ### PROFILE + # PROFILE plt.sca(self.ax_profile) lps = self.line_profile.list_points xa, ya, xb, yb = lps[0].x, lps[0].y, lps[1].x, lps[1].y coor_pix = list( np.asarray( - (np.array([xa, ya, xb, yb]) * self.pp.pp_params.pymses.zoom + 0.5) + (np.array([xa, ya, xb, yb]) * self.pp.params.pymses.zoom + 0.5) * self.shape[0], dtype=int, ) ) xpp, ypp = line(*coor_pix) rho_prof = self.rho_map[ypp, xpp] - xp = ((xpp / float(self.shape[0])) - 0.5) / self.pp.pp_params.pymses.zoom - yp = ((ypp / float(self.shape[1])) - 0.5) / self.pp.pp_params.pymses.zoom + xp = ((xpp / float(self.shape[0])) - 0.5) / self.pp.params.pymses.zoom + yp = ((ypp / float(self.shape[1])) - 0.5) / self.pp.params.pymses.zoom print(xp, yp, xpp, ypp) x = np.sqrt(np.abs((xp - xa) * (xb - xa) + (yp - ya) * (yb - ya))) x = x - float(x[0]) @@ -482,27 +448,27 @@ class InteractiveGUI: self.fit_prof_vmax = vmax self.draw() - def __init__(self, num, path="./", pp_params=None, datamap_key="fluct_coldens_z"): + def __init__(self, num, path="./", params=None, datamap_key="fluct_coldens_z"): """ Interactive plotting """ - if pp_params is None: - pp_params = default_params() - pp_params.input.nml_filename = "disk.nml" - pp_params.out.interactive = True - pp_params.pymses.map_size = 4096 - pp_params.pymses.zoom = 4 + if params is None: + params = default_params() + params.input.nml_filename = "disk.nml" + params.out.interactive = True + params.pymses.map_size = 4096 + params.pymses.zoom = 4 - pp_params.pymses.variables = ["rho", "vel", "P"] + params.pymses.variables = ["rho", "vel", "P"] - pp_params.disk.enable = True - pp_params.disk.nb_bin = 200 - pp_params.pdf.nb_bin = 100 + params.disk.enable = True + params.disk.nb_bin = 200 + params.pdf.nb_bin = 100 self.fig = plt.figure(figsize=(10, 8)) - self.pp = PostProcessor(path, num, pp_params=pp_params, tag="interactive") + self.pp = SnapshotProcessor(path, num, params=params, tag="interactive") self.pp.pdf_coldens("z") self.pp.pdf_rho("z") self.pp.pdf_T("z") @@ -523,8 +489,8 @@ class InteractiveGUI: self.rr = np.sqrt(self.xx ** 2 + self.yy ** 2) self.rhocs = np.column_stack((self.frho_map.flatten(), self.fcs_map.flatten())) - self.rmin = self.pp.pp_params.disk.rmin_pdf / 2.0 - self.rmax = self.pp.pp_params.disk.rmax_pdf / 2.0 + self.rmin = self.pp.params.disk.rmin_pdf / 2.0 + self.rmax = self.pp.params.disk.rmax_pdf / 2.0 self.mask = (self.rr >= self.rmin) & (self.rr <= self.rmax) self.mask_flat = self.mask.flatten() @@ -553,11 +519,11 @@ class InteractiveGUI: ax_rmax = plt.axes([0.05, 0.07, 0.15, 0.02]) self.srmax = Slider( - ax_rmax, r"$r_{max}$", 0, 0.7 / pp_params.pymses.zoom, valinit=self.rmax + ax_rmax, r"$r_{max}$", 0, 0.7 / params.pymses.zoom, valinit=self.rmax ) ax_rmin = plt.axes([0.05, 0.03, 0.15, 0.02]) self.srmin = Slider( - ax_rmin, r"$r_{min}$", 0, 0.7 / pp_params.pymses.zoom, valinit=self.rmin + ax_rmin, r"$r_{min}$", 0, 0.7 / params.pymses.zoom, valinit=self.rmin ) ax_lasso = plt.axes([0.6, 0.07, 0.19, 0.02]) ax_poly = plt.axes([0.8, 0.07, 0.19, 0.02]) diff --git a/ism.py b/ism.py index f513474..bfa6c31 100644 --- a/ism.py +++ b/ism.py @@ -3,7 +3,7 @@ import numpy as np import pandas as pd from plotter import U -import postprocessor +import snapshotprocessor mp = 1.4 * 1.66 * 10 ** (-24) * U.g z0 = 150 * U.pc @@ -53,7 +53,7 @@ def get_sinks(pp): def analyze_box(pp): - pp.cells["mass"] = postprocessor.mass_func(pp.cells) + pp.cells["mass"] = snapshotprocessor.mass_func(pp.cells) pp.coldens("z") coldens = pp.get_value("/maps/coldens_z", unit=U.coldens) sinks = get_sinks(pp) diff --git a/pp_params.py b/params.py similarity index 93% rename from pp_params.py rename to params.py index 6be735d..0270532 100644 --- a/pp_params.py +++ b/params.py @@ -33,4 +33,4 @@ def load_params(filename): def default_params(): - return load_params(_dir_path + "/pp_params.yml") + return load_params(_dir_path + "/params.yml") diff --git a/pp_params.yml b/params.yml similarity index 100% rename from pp_params.yml rename to params.yml diff --git a/pipeline.py b/pipeline.py index 4dd26a1..54d608d 100644 --- a/pipeline.py +++ b/pipeline.py @@ -7,11 +7,11 @@ import os import numpy as np from plotter import Plotter, plt -from postprocessor import PostProcessor, get_time -from comparator import Comparator -from pp_params import default_params, load_params +from snapshotprocessor import SnapshotProcessor, get_time +from studyprocessor import StudyProcessor +from params import default_params, load_params -fake_pp = PostProcessor() +fake_pp = SnapshotProcessor() parser = argparse.ArgumentParser() @@ -192,23 +192,23 @@ storage_in = args.input_path storage_out = args.output_path if args.config is None: - pp_params = default_params() + params = default_params() else: - pp_params = load_params(args.config) + 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 +params.out.zoom = args.zoom +params.out.tag = args.tag +params.out.map_size = args.map_size +params.out.interactive = args.interactive -pp_params.pymses.fft = args.fft +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 +params.disk.on = args.disk +params.disk.binning = args.binning +params.disk.nb_bin = args.nb_bin -pp_params.pdf.nb_bin = args.pdf_nb_bin +params.pdf.nb_bin = args.pdf_nb_bin # extension for out files plt.style.use("seaborn-deep") @@ -265,7 +265,7 @@ for run in runs: while not success: try: if len(args.process) > 0: - pp = PostProcessor(run, num, pp_params=pp_params) + pp = SnapshotProcessor(run, num, params=params) pp.process( args.process, args.process_args, @@ -299,7 +299,7 @@ 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 = Plotter(path_in, runs, nums_success, path_out=path_out, params=params) pl.process( args.plot, args.plot_args, @@ -308,7 +308,7 @@ if len(args.plot) > 0: ) if len(args.compare) > 0: - cc = Comparator(path_in, runs, nums_success, path_out=path_out, pp_params=pp_params) + cc = StudyProcessor(path_in, runs, nums_success, path_out=path_out, params=params) cc.process( args.compare, args.compare_args, diff --git a/plotter.py b/plotter.py index 566b7ad..3f5adeb 100644 --- a/plotter.py +++ b/plotter.py @@ -30,7 +30,7 @@ import matplotlib.pyplot as plt import pspec_read from baseprocessor import Rule, BaseProcessor from aggregator import Aggregator -from comparator import Comparator +from studyprocessor import StudyProcessor from run_selector import RunSelector from units import U, unit_str, convert_exp @@ -121,7 +121,7 @@ class Plotter(Aggregator, BaseProcessor): in_runs=None, in_nums=None, path_out=None, - pp_params=None, + params=None, selector=None, tag=None, unit_time=U.year, @@ -139,19 +139,19 @@ class Plotter(Aggregator, BaseProcessor): in_nums : list or dict of the outputs numbers to consider (ex [3, 5] or {'run1' : [3, 5], 'run2' : [4, 6]) path_out : Path where the plot will be saved. By default set to `path` - pp_params : Parameters for postprocessing. See pp_params module. + params : Parameters for postprocessing. See params module. selector : Existing instance of RunSelector, that selects runs and outputs. If set, in_runs and in_nums will be ignored tag : string to add in the output and data files. kwargs : Keyword arguments for RunSelector. """ - super(Plotter, self).__init__(path, path_out, pp_params, tag) + super(Plotter, self).__init__(path, path_out, params, tag) # Select runs if selector is None: self.selector = RunSelector( - path, in_runs, in_nums, self.pp_params.input.nml_filename, **kwargs + path, in_runs, in_nums, self.params.input.nml_filename, **kwargs ) else: self.selector = selector @@ -161,22 +161,22 @@ class Plotter(Aggregator, BaseProcessor): self.runs = self.selector.runs self.nums = self.selector.nums - # Get comparator object - self.comp = Comparator( + # Get studyprocessor object + self.study = StudyProcessor( path, self.runs, self.nums, path_out, - self.pp_params, + self.params, unit_time=unit_time, selector=self.selector, ) # Get postprocesor objets for each run - self.pp = self.comp.pp + self.snaps = self.study.snaps # Define log prefix - self.log_id = "[plot {}] ".format(self.pp_params.out.tag) + self.log_id = "[plot {}] ".format(self.params.out.tag) # Define rules self.def_rules() @@ -188,12 +188,12 @@ class Plotter(Aggregator, BaseProcessor): def gen_simus(self): self.simulations = {} - simu_fmt = self.pp_params.astrophysix.simu_fmt - descr_fmt = self.pp_params.astrophysix.descr_fmt - tag = self.pp_params.out.tag + simu_fmt = self.params.astrophysix.simu_fmt + descr_fmt = self.params.astrophysix.descr_fmt + tag = self.params.out.tag for run in self.runs: - pp = self.pp[run][self.nums[run][0]] - nml = self.comp.namelist[run] + pp = self.snaps[run][self.nums[run][0]] + nml = self.study.namelist[run] name = simu_fmt.format(run=run, tag=tag, nml=nml) exec_time = str(datetime.datetime.fromtimestamp(os.stat(pp.path).st_ctime)) exec_time = exec_time.split(".")[0] @@ -210,7 +210,7 @@ class Plotter(Aggregator, BaseProcessor): for param in ramses.input_parameters: value = None try: - value = self.comp.get_nml(param.key, run) + value = self.study.get_nml(param.key, run) except KeyError as e: self._log("key {} not found".format(e), "WARNING") @@ -237,8 +237,8 @@ class Plotter(Aggregator, BaseProcessor): """ Check if the dependency belongs to the plotter object or to another one (comp, pp, ..) """ - if dep in self.comp.rules: - result = self.comp.process( + if dep in self.study.rules: + result = self.study.process( dep, dep_arg, overwrite, self.overwrite_dep, select ) if result is not None: @@ -251,7 +251,7 @@ class Plotter(Aggregator, BaseProcessor): Returns true if the plot needs to be redone """ return ( - self.pp_params.out.interactive + self.params.out.interactive or overwrite or not os.path.exists(plot_filename) ) @@ -288,7 +288,7 @@ class Plotter(Aggregator, BaseProcessor): name_full = name # get filetype of the output - filetype = filetype_from_ext[self.pp_params.out.ext] + filetype = filetype_from_ext[self.params.out.ext] # Select runs and nums if select is not None: @@ -330,14 +330,14 @@ class Plotter(Aggregator, BaseProcessor): # Find plot save if from_cells or rule.kind == "cells": - if not os.path.exists(self.pp[run][num].cells_filename): + if not os.exists(self.pp[run][num].cells_filename): self.pp[run][num].load_cells() self.pp[run][num].unload_cells() save = tables.open_file(self.pp[run][num].cells_filename) elif rule.kind == "snapshot": - save = tables.open_file(self.pp[run][num].filename) + save = tables.open_file(self.snaps[run][num].filename) else: - save = tables.open_file(self.comp.filename, "r") + save = tables.open_file(self.study.filename, "r") # Call plot routine try: @@ -362,7 +362,7 @@ class Plotter(Aggregator, BaseProcessor): if plot_info is not None: df.plot_info = plot_info if num is not None: - snap = self.pp[run][num].snapshot + snap = self.snaps[run][num].snapshot if overwrite and df.name in snap.datafiles: del snap.datafiles[df.name] @@ -385,10 +385,10 @@ class Plotter(Aggregator, BaseProcessor): if self._needs_computation(overwrite, plot_filename): plot_info = rule.plot(save, arg, **kwargs) - if not self.pp_params.out.interactive and close: + if not self.params.out.interactive and close: plt.tight_layout(pad=1) - if self.pp_params.out.save: + if self.params.out.save: plt.savefig(plot_filename) self._log("{} plotted".format(plot_filename), "SUCCESS") else: @@ -396,7 +396,7 @@ class Plotter(Aggregator, BaseProcessor): "{} plotted".format(os.path.basename(plot_filename)), "SUCCESS" ) - if not self.pp_params.out.interactive and close: + if not self.params.out.interactive and close: plt.close() return plot_info else: @@ -406,10 +406,10 @@ class Plotter(Aggregator, BaseProcessor): """ Determine a filename based on rule name, run, output and parameters """ - tag_name = self.pp_params.out.tag + tag_name = self.params.out.tag - if fmt is None and self.pp_params.out.fmt == "": - if not self.pp_params.out.tag == "": + if fmt is None and self.params.out.fmt == "": + if not self.params.out.tag == "": tag_name = "_" + tag_name if run is not None and num is not None: @@ -419,11 +419,11 @@ class Plotter(Aggregator, BaseProcessor): else: fmt = "{out}/{name}{tag}{ext}" elif fmt is None: - fmt = self.pp_params.out.fmt + fmt = self.params.out.fmt nml = None if run is not None: - nml = self.comp.namelist[run] + nml = self.study.namelist[run] return fmt.format( run=run, @@ -432,7 +432,7 @@ class Plotter(Aggregator, BaseProcessor): num=num, nml=nml, out=self.path_out, - ext=self.pp_params.out.ext, + ext=self.params.out.ext, ) def get_label_run(self, run, label=None, nml_key=None, time=None): @@ -514,11 +514,9 @@ class Plotter(Aggregator, BaseProcessor): title = self.get_label_run(run, title, nml_key) if put_time: - time = self.save.root._v_attrs.time * self.comp.info["unit_time"] + time = self.save.root._v_attrs.time * self.study.info["unit_time"] u_str = unit_str(unit_time, format="{unit}") - time_str = self.pp_params.plot.time_fmt.format( - time.express(unit_time), u_str - ) + time_str = self.params.plot.time_fmt.format(time.express(unit_time), u_str) if len(title) > 0: title = title + " | " + time_str else: @@ -593,7 +591,7 @@ class Plotter(Aggregator, BaseProcessor): dmap, extent=im_extent, origin="lower", norm=norm, cmap=cmap, **kwargs ) - plt.locator_params(axis="both", nbins=self.pp_params.plot.ntick) + plt.locator_params(axis="both", nbins=self.params.plot.ntick) if xlabel is None: xlabel = self._ax_title[ax_h] @@ -781,7 +779,7 @@ class Plotter(Aggregator, BaseProcessor): label, unit_old, unit = self._ax_label_unit(dmap_vh_node, "", unit, unit_coeff) - vel_red = self.pp_params.plot.vel_red + vel_red = self.params.plot.vel_red # take only a subset of velocities map_vh_red = dmap_vh[::vel_red, ::vel_red] * unit_old.express(unit) @@ -835,7 +833,7 @@ class Plotter(Aggregator, BaseProcessor): # TODO : redo this with im_extent - vel_red = self.pp_params.plot.vel_red + vel_red = self.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 @@ -923,7 +921,7 @@ class Plotter(Aggregator, BaseProcessor): if nml_color is None: color = colors[run] else: - nml = self.comp.get_nml(nml_color, run) + nml = self.study.get_nml(nml_color, run) try: color = colors[nml] except TypeError: @@ -1043,8 +1041,8 @@ class Plotter(Aggregator, BaseProcessor): # If relevent, get time if put_time: - time = self.save.root._v_attrs.time * self.comp.info["unit_time"] - time_str = self.pp_params.plot.time_fmt.format( + time = self.save.root._v_attrs.time * self.study.info["unit_time"] + time_str = self.params.plot.time_fmt.format( time.express(unit_time), unit_time.latex.replace("text", "math") ) time_str = f"${time_str}$" @@ -1125,11 +1123,11 @@ class Plotter(Aggregator, BaseProcessor): color = colors[run] elif nml_color == "time": time = ( - self.save.root._v_attrs.time * self.comp.info["unit_time"] + self.save.root._v_attrs.time * self.study.info["unit_time"] ).express(unit_time) color = colors(time) else: - nml = self.comp.get_nml(nml_color, run) + nml = self.study.get_nml(nml_color, run) try: color = colors[nml] except TypeError: @@ -1254,7 +1252,7 @@ class Plotter(Aggregator, BaseProcessor): ssfr_sun = 2.5e-9 ssfr_ken = ssfr_sun * n0 ** 1.4 - coeff = ssfr_ken * 1e6 * (self.comp.info["unit_length"].express(U.pc)) ** 2 + coeff = ssfr_ken * 1e6 * (self.study.info["unit_length"].express(U.pc)) ** 2 for i in np.arange(tmin, max(tmax, tmin + ymax / coeff), step): t = np.linspace(0, tmax, 1000) plt.plot(t + i, t * coeff, ls="--", lw=0.9, color="grey") @@ -1661,7 +1659,7 @@ class Plotter(Aggregator, BaseProcessor): ] # Generic rules directly from Ramses fields - for field in self.pp_params.pymses.variables: + for field in self.params.pymses.variables: def generic_rule(name): diff --git a/postprocessor.py b/snapshotprocessor.py similarity index 93% rename from postprocessor.py rename to snapshotprocessor.py index 3472017..dfe24f9 100644 --- a/postprocessor.py +++ b/snapshotprocessor.py @@ -42,8 +42,7 @@ from baseprocessor import ( oct_vect_getter, ) -from run_selector import NamelistRecursive -import f90nml +from run_selector import RunSelector # Getters @@ -255,7 +254,8 @@ def degrade_map(dmap, nnew, integrate=False): return dmap_new -class PostProcessor(HDF5Container): +class SnapshotProcessor(HDF5Container): + """ This class enable to compute and save derived quantities from the raw output """ @@ -286,19 +286,29 @@ class PostProcessor(HDF5Container): path=None, num=None, path_out=None, - pp_params=None, + params=None, tag=None, + selector=None, unit_time=U.year, ): """ Creates the basic structures needed for the outputs + + Parameters + ---------- + path : str, path of where the outputs are located + path_out : str, where to store postprocessings results + params : str or params instance + tag : str, distinguish this postprocessing set + selector : RunSelector instance + unit_time : Unit instance, used for astrophysix """ - super(PostProcessor, self).__init__(path, path_out, pp_params, tag) + super(SnapshotProcessor, self).__init__(path, path_out, params, tag) # Open outfile - if not self.pp_params.out.tag == "": - tag_name = self.pp_params.out.tag + "_" + if not self.params.out.tag == "": + tag_name = self.params.out.tag + "_" else: tag_name = "" @@ -325,113 +335,34 @@ class PostProcessor(HDF5Container): if not os.path.exists(self.path_out): os.makedirs(self.path_out) - self.open() - - # Ramses Output self.path = path self.run = os.path.basename(path) self.num = num - self._ro = pymses.RamsesOutput( - path, - num, - order=self.pp_params.pymses.order, - verbose=self.pp_params.pymses.verbose, - ) - self._amr = self._ro.amr_source(self.pp_params.pymses.variables) - self._part = self._ro.particle_source(self.pp_params.pymses.part_variables) - if self.pp_params.pymses.filter: - self.min_coords = np.array(self.pp_params.pymses.min_coords) - self.max_coords = np.array(self.pp_params.pymses.max_coords) - box = reg.Box((self.min_coords, self.max_coords)) - - amr_filt = RegionFilter(box, self._amr) - self._amr = amr_filt - part_filt = RegionFilter(box, self._part) - self._part = part_filt - - self.info = self._ro.info.copy() - - # Density operator - self._rho_op = ScalarOperator( - lambda dset: dset["rho"], self._ro.info["unit_density"] - ) - - # Density ray tracer - if self.pp_params.pymses.fft: - self._rt = splatting.SplatterProcessor( - self._amr, self._ro.info, self._rho_op + # Create selector object + if selector is None: + selector = RunSelector( + os.path.dirname(path), + self.run, + self.num, + self.params.input.nml_filename, ) - else: - self._rt = raytracing.RayTracer(self._amr, self._ro.info, self._rho_op) - if not self.pp_params.pymses.multiprocessing: - self._rt.disable_multiprocessing() - - # Set the extent of the image - self._radius = 0.5 / self.pp_params.pymses.zoom + self.info = selector.info[self.run][self.num] + self.namelist = selector.namelist[self.run][self.num] self.lbox = self.info["boxlen"] - if self.pp_params.pymses.filter: - center = (self.max_coords + self.min_coords) / 2.0 - im_extent = np.array( - [ - self.min_coords[0], - self.max_coords[0], - self.min_coords[1], - self.max_coords[1], - ] - ) - distance = (self.max_coords[2] - self.min_coords[2]) / 2.0 - else: - center = self.pp_params.pymses.center - im_extent = np.array( - [ - (-self._radius + center[0]), - (self._radius + center[0]), - (-self._radius + center[1]), - (self._radius + center[1]), - ] - ) - distance = self._radius - # Get time - self.time = self._ro.info["time"] - - # Get namelist - self.namelist = None - path_nml = path + "/" + self.pp_params.input.nml_filename - self.namelist = NamelistRecursive(f90nml.read(path_nml)) + self.time = self.info["time"] # Set post processing attributes + self.open() 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.unit_length = self.info["unit_length"] self.save.root._v_attrs.time = self.time - - if "/maps" not 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 = {} - for ax_los in self._ax_nb: # los = line of sight - ax_v = self._axes_v[ax_los] - - self._cam[ax_los] = Camera( - center=self.pp_params.pymses.center, - line_of_sight_axis=ax_los, - region_size=[im_extent[1] - im_extent[0], im_extent[3] - im_extent[2]], - distance=distance, - far_cut_depth=distance, - up_vector=ax_v, - map_max_size=self.pp_params.pymses.map_size, - ) - self.close() self.log_id = "[{}, {}] ".format(self.run, self.num) @@ -451,8 +382,97 @@ class PostProcessor(HDF5Container): data_reference="OUTPUT_{}".format(self.num), ) + try: + self.init_pymses() + except FileNotFoundError: + self._log("Pymses not initialized", "WARNING") + self.def_rules() + def init_pymses(self): + self._ro = pymses.RamsesOutput( + self.path, + self.num, + order=self.params.pymses.order, + verbose=self.params.pymses.verbose, + ) + self._amr = self._ro.amr_source(self.params.pymses.variables) + self._part = self._ro.particle_source(self.params.pymses.part_variables) + + if self.params.pymses.filter: + self.min_coords = np.array(self.params.pymses.min_coords) + self.max_coords = np.array(self.params.pymses.max_coords) + box = reg.Box((self.min_coords, self.max_coords)) + + amr_filt = RegionFilter(box, self._amr) + self._amr = amr_filt + part_filt = RegionFilter(box, self._part) + self._part = part_filt + + # Density operator + self._rho_op = ScalarOperator( + lambda dset: dset["rho"], self._ro.info["unit_density"] + ) + + # Density ray tracer + if self.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) + + if not self.params.pymses.multiprocessing: + self._rt.disable_multiprocessing() + + # Set the extent of the image + self._radius = 0.5 / self.params.pymses.zoom + if self.params.pymses.filter: + center = (self.max_coords + self.min_coords) / 2.0 + im_extent = np.array( + [ + self.min_coords[0], + self.max_coords[0], + self.min_coords[1], + self.max_coords[1], + ] + ) + distance = (self.max_coords[2] - self.min_coords[2]) / 2.0 + else: + center = self.params.pymses.center + im_extent = np.array( + [ + (-self._radius + center[0]), + (self._radius + center[0]), + (-self._radius + center[1]), + (self._radius + center[1]), + ] + ) + distance = self._radius + + # Initialize cameras + self._cam = {} + for ax_los in self._ax_nb: # los = line of sight + ax_v = self._axes_v[ax_los] + + self._cam[ax_los] = Camera( + center=self.params.pymses.center, + line_of_sight_axis=ax_los, + region_size=[im_extent[1] - im_extent[0], im_extent[3] - im_extent[2]], + distance=distance, + far_cut_depth=distance, + up_vector=ax_v, + map_max_size=self.params.pymses.map_size, + ) + + self.open() + if "/maps" not 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 + self.close() + def load_data(self, points_src, filename, save, keys=None): """ Load data from the source file in the memory. @@ -496,7 +516,7 @@ class PostProcessor(HDF5Container): self.parts = self.load_data( self._part, self.parts_filename, - self.pp_params.process.save_parts, + self.params.process.save_parts, keys=keys, ) self.parts_loaded = True @@ -521,7 +541,7 @@ class PostProcessor(HDF5Container): self.cells = self.load_data( cells_src, self.cells_filename, - self.pp_params.process.save_cells, + self.params.process.save_cells, keys=keys, ) self.cells_loaded = True @@ -551,7 +571,7 @@ class PostProcessor(HDF5Container): pos = dset.points except AttributeError: pos = dset["pos"] - pos = pos - np.array(self.pp_params.disk.center) + pos = pos - np.array(self.params.disk.center) return pos def getter_vect_r(self, dset, name_vect): @@ -576,7 +596,7 @@ class PostProcessor(HDF5Container): Returns the position in normalized units centered on the position of the star """ pos = dset.get_cell_centers() - pos = pos - np.array(self.pp_params.disk.center) + pos = pos - np.array(self.params.disk.center) return pos def oct_getter_vect_r(self, dset, name_vect): @@ -665,12 +685,12 @@ class PostProcessor(HDF5Container): op = FractionOperator(num, denum, unit) - if self.pp_params.pymses.fft: + if self.params.pymses.fft: rt = splatting.SplatterProcessor(self._amr, self._ro.info, op) else: rt = raytracing.RayTracer(self._amr, self._ro.info, op) - if not self.pp_params.pymses.multiprocessing: + if not self.params.pymses.multiprocessing: rt.disable_multiprocessing() datamap = rt.process(self._cam[ax_los], surf_qty=surf_qty) @@ -705,7 +725,7 @@ class PostProcessor(HDF5Container): else: df["value"] = value - if self.pp_params.process.unload_cells: + if self.params.process.unload_cells: self.unload_cells() df.sort_values("axis", inplace=True) @@ -726,7 +746,7 @@ class PostProcessor(HDF5Container): else: data = np.sum(value, axis=0) - if self.pp_params.process.unload_cells: + if self.params.process.unload_cells: self.unload_cells() return data @@ -744,7 +764,7 @@ class PostProcessor(HDF5Container): # Transpose (.T) is for vectorial values data = np.sum((weight * value.T).T, axis=0) / np.sum(weight) - if self.pp_params.process.unload_cells: + if self.params.process.unload_cells: self.unload_cells() return data @@ -754,7 +774,7 @@ class PostProcessor(HDF5Container): if logbins: data = np.log10(data) weights = weight_func(self.cells) - if self.pp_params.process.unload_cells: + if self.params.process.unload_cells: self.unload_cells() values, edges = np.histogram(data, bins, weights=weights) @@ -771,7 +791,7 @@ class PostProcessor(HDF5Container): centers, B_mean = mean_by_bins(rho, B, bins, logbins) - if self.pp_params.process.unload_cells: + if self.params.process.unload_cells: self.unload_cells() return ({"rho": centers, "B": B_mean}, {"logbins": logbins}) @@ -840,7 +860,7 @@ class PostProcessor(HDF5Container): else: centers = 0.5 * (rho_bins[1:] + rho_bins[:-1]) - if self.pp_params.process.unload_cells: + if self.params.process.unload_cells: self.unload_cells() return ({"rho": centers, "Ek_Eb_rho": ek_eb}, {"logbins": logbins}) @@ -928,10 +948,10 @@ class PostProcessor(HDF5Container): return dmap_P / dmap_rho def _levels(self, ax_los): - self._amr.set_read_levelmax(self.pp_params.pymses.levelmax) + self._amr.set_read_levelmax(self.params.pymses.levelmax) level_op = MaxLevelOperator() rt_level = raytracing.RayTracer(self._amr, self._ro.info, level_op) - if not self.pp_params.pymses.multiprocessing: + if not self.params.pymses.multiprocessing: rt_level.disable_multiprocessing() datamap = rt_level.process(self._cam[ax_los], surf_qty=True) return datamap.map.T @@ -978,7 +998,7 @@ class PostProcessor(HDF5Container): # Ray tracer for the angular speed rt_omega = raytracing.RayTracer(self._amr, self._ro.info, omega_op) - if not self.pp_params.pymses.multiprocessing: + if not self.params.pymses.multiprocessing: rt_omega.disable_multiprocessing() dmap_omega = rt_omega.process(self._cam[ax_los]).map.T @@ -1020,7 +1040,7 @@ class PostProcessor(HDF5Container): # Get coordinates im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * space_coeff - center = np.array(self.pp_params.disk.center) * space_coeff + center = np.array(self.params.disk.center) * space_coeff # Physical size of cells dx = (im_extent[1] - im_extent[0]) / map_size @@ -1046,7 +1066,7 @@ class PostProcessor(HDF5Container): """ Computes radial bins (for disk) """ - center = np.array(self.pp_params.disk.center) * self.lbox + center = np.array(self.params.disk.center) * self.lbox im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * self.lbox # radius of the corner of the box plus a margin @@ -1055,20 +1075,20 @@ class PostProcessor(HDF5Container): + 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 + bin_in = self.params.disk.bin_in + bin_out = self.params.disk.bin_out + nb_bin = self.params.disk.nb_bin # radial bins - if self.pp_params.disk.binning in ["log", "logarithmic"]: + if self.params.disk.binning in ["log", "logarithmic"]: lrad_in = np.log10(bin_in) lrad_ext = np.log10(bin_out) rad_bins = np.logspace(lrad_in, lrad_ext, num=nb_bin) - elif self.pp_params.disk.binning in ["lin", "linear"]: + elif self.params.disk.binning in ["lin", "linear"]: rad_bins = np.linspace(bin_in, bin_out, num=nb_bin) else: raise RuntimeWarning( - f"Invalid parameter {self.pp_params.disk.binning} for disk binning method, using linear binning" + f"Invalid parameter {self.params.disk.binning} for disk binning method, using linear binning" ) rad_bins = np.linspace(bin_in, bin_out, num=nb_bin) @@ -1086,8 +1106,8 @@ class PostProcessor(HDF5Container): Computes the radius from the center """ im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * self.lbox - map_size = self.pp_params.pymses.map_size - center = np.array(self.pp_params.disk.center) * self.lbox + map_size = self.params.pymses.map_size + center = np.array(self.params.disk.center) * self.lbox # Physical size of cells dx = (im_extent[1] - im_extent[0]) / map_size @@ -1220,15 +1240,15 @@ class PostProcessor(HDF5Container): rr = self.get_value("/maps/rr_" + ax_los) mask_pdf = ( - (rr > self.pp_params.disk.rmin_pdf) - & (rr < self.pp_params.disk.rmax_pdf) + (rr > self.params.disk.rmin_pdf) + & (rr < self.params.disk.rmax_pdf) & (fluct_map > 0) ) values, edges = np.histogram( np.log10(fluct_map[mask_pdf].flatten()), - self.pp_params.pdf.nb_bin, - range=self.pp_params.pdf.range, + self.params.pdf.nb_bin, + range=self.params.pdf.range, density=True, ) centers = 0.5 * (edges[1:] + edges[:-1]) @@ -1238,9 +1258,9 @@ class PostProcessor(HDF5Container): pdf = self.get_value("/hist/pdf_" + name + "_" + ax_los) values, centers = pdf mask_fit = ( - (centers > self.pp_params.pdf.xmin_fit) - & (centers < self.pp_params.pdf.xmax_fit) - & (values > self.pp_params.pdf.fit_cut * np.max(values)) + (centers > self.params.pdf.xmin_fit) + & (centers < self.params.pdf.xmax_fit) + & (values > self.params.pdf.fit_cut * np.max(values)) ) (slope, origin, correlation, _, stderr) = linregress( centers[mask_fit], np.log10(values[mask_fit]) @@ -1315,7 +1335,7 @@ class PostProcessor(HDF5Container): def getter_alpha_grav(dset): r2 = np.sum((self.lbox * self.oct_getter_pos_disk(dset)) ** 2, axis=2) e2 = (1.0 / 256.0) ** 2 - gstar = -self.G * self.pp_params.disk.mass_star / (e2 + r2) + gstar = -self.G * self.params.disk.mass_star / (e2 + r2) gr = self.oct_getter_vect_r(dset, "g") - gstar gphi = self.oct_getter_vect_phi(dset, "g") return gr * gphi / (4 * np.pi * self.G) @@ -1388,14 +1408,14 @@ class PostProcessor(HDF5Container): def _filaments(self): - datamap_name = self.pp_params.filaments.datamap - verbose = self.pp_params.filaments.verbose - rmin_frac = self.pp_params.filaments.rmin - rmax_frac = self.pp_params.filaments.rmax - size_thresh = self.pp_params.filaments.size_thresh - skel_thresh = self.pp_params.filaments.skel_thresh - branch_thresh = self.pp_params.filaments.branch_thresh - glob_thresh = self.pp_params.filaments.glob_thresh + datamap_name = self.params.filaments.datamap + verbose = self.params.filaments.verbose + rmin_frac = self.params.filaments.rmin + rmax_frac = self.params.filaments.rmax + size_thresh = self.params.filaments.size_thresh + skel_thresh = self.params.filaments.skel_thresh + branch_thresh = self.params.filaments.branch_thresh + glob_thresh = self.params.filaments.glob_thresh datamap = self.get_value("/maps/" + datamap_name + "_z") shape = datamap.shape @@ -1451,7 +1471,7 @@ class PostProcessor(HDF5Container): Compute forces within a filament (for disks), within the slice at z=0 """ - GM = self.G * self.pp_params.disk.mass_star # Mass parameter + GM = self.G * self.params.disk.mass_star # Mass parameter # Find center of filaments i_center, j_center = self._filaments_center() @@ -1465,8 +1485,8 @@ class PostProcessor(HDF5Container): # Get coordinates im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * self.lbox - map_size = self.pp_params.pymses.map_size - center = np.array(self.pp_params.disk.center) * self.lbox + map_size = self.params.pymses.map_size + center = np.array(self.params.disk.center) * self.lbox # Physical size of cells dx = (im_extent[1] - im_extent[0]) / map_size @@ -1661,7 +1681,7 @@ class PostProcessor(HDF5Container): self._filaments, "Filaments", "/datasets", - dependencies={self.pp_params.filaments.datamap: "z"}, + dependencies={self.params.filaments.datamap: "z"}, ), "filaments_forces": Rule( self, @@ -1837,7 +1857,7 @@ class PostProcessor(HDF5Container): ] # Generic rules directly from Ramses fields - for field in self.pp_params.pymses.variables: + for field in self.params.pymses.variables: def generic_rule(name, getter, unit, oct_getter=None): if oct_getter is None: @@ -2002,7 +2022,7 @@ class PostProcessor(HDF5Container): self._gen_rule_transform("fluct_coldens", np.nanmax, "max", group="/globals") - super(PostProcessor, self).def_rules() + super(SnapshotProcessor, self).def_rules() def get_time(path, num): diff --git a/comparator.py b/studyprocessor.py similarity index 94% rename from comparator.py rename to studyprocessor.py index 05a377c..1f67623 100644 --- a/comparator.py +++ b/studyprocessor.py @@ -8,15 +8,15 @@ from scipy.stats import linregress from baseprocessor import Rule, HDF5Container from aggregator import Aggregator -from postprocessor import PostProcessor +from snapshotprocessor import SnapshotProcessor from run_selector import RunSelector -from pp_params import default_params +from params import default_params from units import U -class Comparator(Aggregator, HDF5Container): +class StudyProcessor(Aggregator, HDF5Container): """ - Do comparaison between outputs and runs + This object is linked to several ramses simulation of a same study """ def __init__( @@ -25,7 +25,7 @@ class Comparator(Aggregator, HDF5Container): in_runs, in_nums, path_out=None, - pp_params=default_params(), + params=default_params(), selector=None, tag=None, unit_time=U.year, @@ -35,11 +35,11 @@ class Comparator(Aggregator, HDF5Container): Creates the basic structures needed for the outputs """ - super(Comparator, self).__init__(path, path_out, pp_params, tag) + super(StudyProcessor, self).__init__(path, path_out, params, tag) # Open outfile - if not self.pp_params.out.tag == "": - tag_name = "_" + self.pp_params.out.tag + if not self.params.out.tag == "": + tag_name = "_" + self.params.out.tag else: tag_name = "" @@ -48,7 +48,7 @@ class Comparator(Aggregator, HDF5Container): # Select runs if selector is None: selector = RunSelector( - path, in_runs, in_nums, self.pp_params.input.nml_filename, **kwargs + path, in_runs, in_nums, self.params.input.nml_filename, **kwargs ) # Save infos @@ -57,30 +57,28 @@ class Comparator(Aggregator, HDF5Container): self.nums = selector.nums # Get postprocesor objets for each run and infos on them - self.pp = {} + self.snaps = {} self.info = {} for run in self.runs: path_run = path + "/" + run path_out_run = path_out + "/" + run - self.pp[run] = {} + self.snaps[run] = {} for num in self.nums[run]: - self.pp[run][num] = PostProcessor( + self.snaps[run][num] = SnapshotProcessor( path_run, num, path_out=path_out_run, - pp_params=self.pp_params, + params=self.params, unit_time=unit_time, ) run0 = self.runs[0] - self.info = self.pp[run0][self.nums[run0][0]].info - + self.info = selector.info[run0][self.nums[run0][0]] self.namelist = selector.namelist - # log info - self.log_id = "[comp {}] ".format(self.pp_params.out.tag) + self.log_id = "[comp {}] ".format(self.params.out.tag) # Define rules self.def_rules() @@ -107,7 +105,7 @@ class Comparator(Aggregator, HDF5Container): return missing_nums def _save_data(self, name_full, data, description, unit): - super(Comparator, self)._save_data(name_full, data, description, unit) + super(StudyProcessor, self)._save_data(name_full, data, description, unit) self.save.get_node(name_full)._v_attrs.nums = self.nums def _time_series(self, getter, arg=None): @@ -194,13 +192,13 @@ class Comparator(Aggregator, HDF5Container): } def get_attr(self, attr_name, run, num, node_name="/", arg=None): - pp = self.pp[run][num] + pp = self.snaps[run][num] if arg is not None: node_name = node_name + "_" + str(arg) return pp.get_attribute(node_name, attr_name) - def get_pp_value(self, name, run, num, arg=None): - pp = self.pp[run][num] + def get_snap_value(self, name, run, num, arg=None): + pp = self.snaps[run][num] if arg is not None: name = name + "_" + str(arg) return pp.get_value(name) @@ -208,7 +206,7 @@ class Comparator(Aggregator, HDF5Container): def get_global(self, node_name, run, num, arg=None, unload_cells=False): if arg is not None: node_name = node_name + "_" + str(arg) - pp = self.pp[run][num] + pp = self.snaps[run][num] if unload_cells: pp.unload_cells() value = pp.get_value(node_name) @@ -218,7 +216,7 @@ class Comparator(Aggregator, HDF5Container): return self.namelist[run][nml_key] def get_pdf_slope(self, name, run, num): - pp = self.pp[run][num] + pp = self.snaps[run][num] pp.process(["fit_pdf_" + name], ["z"], overwrite=self.overwrite_dep) slope = pp.get_attribute("/hist/pdf_" + name + "_z", "slope") return slope @@ -274,7 +272,7 @@ class Comparator(Aggregator, HDF5Container): return series def _extract_coarse_step_from_log(self, series, log_filename, run): - rism = self.pp_params.input.ramses_ism + rism = self.params.input.ramses_ism nlines = 2 + int(rism) # Number of useful lines cmd_grep = "grep 'Main step\\|coarse step' {} -A {}".format( log_filename, nlines - 1 @@ -324,7 +322,7 @@ class Comparator(Aggregator, HDF5Container): series["turb_rms"][run].append(np.float(content[i + 1].split(":")[1])) try: turb_energy = np.float(content[i + 2].split(":")[1]) - threshold = self.pp_params.rules.turb_energy_threshold + threshold = self.params.rules.turb_energy_threshold assert threshold < 0 or abs(turb_energy) < threshold series["turb_energy"][run].append(abs(turb_energy)) except (AssertionError, ValueError, IndexError): @@ -347,7 +345,7 @@ class Comparator(Aggregator, HDF5Container): path_run = self.path + "/" + run # Get list of run files - log_files = path_run + "/" + self.pp_params.input.log_prefix + "*" + log_files = path_run + "/" + self.params.input.log_prefix + "*" # Parse files for log_filename in glob.glob(log_files): @@ -372,7 +370,7 @@ class Comparator(Aggregator, HDF5Container): ssfr = {} for run in self.runs: # Surface of the box in pc^2 - info = self.pp[run][self.nums[run][0]].info + info = self.snaps[run][self.nums[run][0]].info surface = (info["unit_length"].express(U.pc)) ** 2 # WARNING : We do not multiply by boxlen since already done in 'unit_length' (pymses) @@ -404,7 +402,7 @@ class Comparator(Aggregator, HDF5Container): ssm = {} for run in self.runs: # Surface of the box in pc^2 - info = self.pp[run][self.nums[run][0]].info + info = self.snaps[run][self.nums[run][0]].info surface = (info["unit_length"].express(U.pc)) ** 2 mass_sink = self.save.get_node( "/series/sinks_from_log/mass_sink/" + run @@ -441,9 +439,9 @@ class Comparator(Aggregator, HDF5Container): for i, run in enumerate(col_pdf["runs"]): values, centers = col_pdf["mean"][i] mask_fit = ( - (centers > self.pp_params.pdf.xmin_fit) - & (centers < self.pp_params.pdf.xmax_fit) - & (values > np.max(values) * self.pp_params.pdf.fit_cut) + (centers > self.params.pdf.xmin_fit) + & (centers < self.params.pdf.xmax_fit) + & (values > np.max(values) * self.params.pdf.fit_cut) ) (slope[i], origin[i], correlation, _, stderr[i]) = linregress( centers[mask_fit], np.log10(values[mask_fit]) @@ -740,4 +738,4 @@ class Comparator(Aggregator, HDF5Container): self._gen_rule_avg("rms_from_log", "turb_rms") self._gen_rule_avg("rms_from_log", "turb_energy") - super(Comparator, self).def_rules() + super(StudyProcessor, self).def_rules()