diff --git a/aggregator.py b/aggregator.py index 02d3fec..aa02f90 100644 --- a/aggregator.py +++ b/aggregator.py @@ -1,5 +1,11 @@ +# coding: utf-8 + +import numpy as np +from functools import partial from mypool import MyPool -from postprocessor import * + +from postprocessor import PostProcessor +from run_selector import RunSelector def _map_rule(rule, arg, overwrite, path, path_out, pp_params, run_num): diff --git a/baseprocessor.py b/baseprocessor.py index bf9a147..4185428 100644 --- a/baseprocessor.py +++ b/baseprocessor.py @@ -1,44 +1,18 @@ # coding: utf-8 import copy -import glob as glob import os -import subprocess -import sys import time -from abc import ABCMeta, abstractmethod +from abc import ABCMeta from functools import partial import numpy as np -import pymses -import tables -from astrophysix.simdm import SimulationStudy -from astrophysix.simdm.experiment import ( - AppliedAlgorithm, - ParameterSetting, - ParameterVisibility, - Simulation, -) -from astrophysix.simdm.results import GenericResult, Snapshot -from numpy.polynomial.polynomial import polyfit -from pymses.analysis import ( - Camera, - FractionOperator, - MaxLevelOperator, - ScalarOperator, - raytracing, - slicing, - splatting, -) -from pymses.filters import CellsToPoints -from pymses.sources.hop.file_formats import * -from pymses.sources.ramses import output -from scipy.stats import linregress -from tables import HDF5ExtError -from ramses_astrophysix import ramses -from run_selector import * -from units import * +import tables + +from tables import HDF5ExtError +from pp_params import default_params, load_params +from units import U class Rule: @@ -51,7 +25,7 @@ class Rule: dependencies=[], is_valid=lambda arg: True, kind="classic", - unit=cst.none, + unit=U.none, ): self.postproc = postproc self.process_fn = process @@ -63,24 +37,11 @@ class Rule: self.kind = kind def process(self, arg, **kwargs): - if not arg is None: + if arg is not None: return self.process_fn(arg, **kwargs) else: return self.process_fn(**kwargs) - def is_valid(self, arg): - # save = self.postproc.save - # valid = True - # for dep in self.dependencies: - # if dep in self.postproc.rules: - # 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 valid and self.is_valid_add(arg) - return self.is_valid_add(arg) - class BaseProcessor: """ @@ -134,7 +95,9 @@ class BaseProcessor: ) else: self._log( - "{} is unknown, allowed rules are {}".format(name, self.rules.keys()), + "{} is unknown, allowed rules are {}".format( + to_process, self.rules.keys() + ), "ERROR", ) @@ -176,13 +139,13 @@ class BaseProcessor: return overwrite def _process_rule(self, name, rule, arg, overwrite=False, **kwargs): - if not arg is None: + if arg is not None: name_full = rule.group + "/" + name + "_" + str(arg) else: name_full = rule.group + "/" + name if rule.is_valid(arg): - if not name_full in self.just_done: + if name_full not in self.just_done: if self._needs_computation(overwrite, name_full): self._log("Processing {}".format(name_full)) data = rule.process(arg, **kwargs) @@ -255,7 +218,7 @@ class HDF5Container(BaseProcessor): ) else: value = node.read() - if not (unit is None or unit_old is None or unit_old == cst.none): + if not (unit is None or unit_old is None or unit_old == U.none): value = value * unit_old.express(unit) finally: if not open_before: @@ -266,7 +229,7 @@ class HDF5Container(BaseProcessor): """ Get real units from info files unit is either: - 1. An instance of cst.Unit (pymses unit class) + 1. An instance of U.Unit (pymses unit class) 2. A string beginning by "unit_", referring to a code unit, available in self.info 3. A dict {unit1 : exp1, unit2: exp2, ...} with unitX as 2. @@ -276,10 +239,10 @@ class HDF5Container(BaseProcessor): and unit the corresponding unit (on one on the above format) Returns: - 1-3. : a cst.Unit instance - 4. : a dict {key: unit, ...} with same key as input and unit being cst.Unit instances + 1-3. : a U.Unit instance + 4. : a dict {key: unit, ...} with same key as input and unit being U.Unit instances """ - if isinstance(unit, cst.Unit): + if isinstance(unit, U.Unit): return unit if isinstance(unit, str) and unit[:5] == "unit_": res = self.info[unit] @@ -287,13 +250,13 @@ class HDF5Container(BaseProcessor): res = res / self.info["boxlen"] return res if list(unit)[0][:5] == "unit_": - new_unit = cst.none + new_unit = U.none for base_unit_str in unit: expo = unit[base_unit_str] base_unit = self._get_units(base_unit_str) new_unit = new_unit * base_unit ** expo return new_unit - if (not data is None) and isinstance(data, dict) and list(unit)[0] in data: + if (data is not None) and isinstance(data, dict) and list(unit)[0] in data: for key in unit: unit[key] = self._get_units(unit[key]) return unit @@ -364,7 +327,7 @@ class HDF5Container(BaseProcessor): ) self.save.get_node(name_full).attrs.unit = unit - if not attrs is None: + if attrs is not None: for key in attrs: key = str(key) self.save.get_node(name_full)._v_attrs[key] = attrs[key] @@ -448,7 +411,7 @@ class HDF5Container(BaseProcessor): fn, name_array_in, name_array_out, - unit_out=cst.none, + unit_out=U.none, description="", recursive=True, ): diff --git a/comparator.py b/comparator.py index ee050ef..e8d5bb9 100644 --- a/comparator.py +++ b/comparator.py @@ -1,6 +1,17 @@ # coding: utf-8 -from aggregator import * +import os +import glob +import numpy as np +from functools import partial +from scipy.stats import linregress + +from baseprocessor import Rule, HDF5Container +from aggregator import Aggregator +from postprocessor import PostProcessor +from run_selector import RunSelector +from pp_params import default_params +from units import U class Comparator(Aggregator, HDF5Container): @@ -17,7 +28,7 @@ class Comparator(Aggregator, HDF5Container): pp_params=default_params(), selector=None, tag=None, - unit_time=cst.year, + unit_time=U.year, **kwargs ): """ @@ -78,14 +89,14 @@ class Comparator(Aggregator, HDF5Container): """ if overwrite or not (name_full in self.save): return True - elif not "nums" in self.save.get_node(name_full)._v_attrs: + elif "nums" not in self.save.get_node(name_full)._v_attrs: return True else: saved_nums = self.save.get_node(name_full)._v_attrs.nums - missing_runs = len([run for run in self.nums if not run in saved_nums]) > 0 + missing_runs = len([run for run in self.nums if run not in saved_nums]) > 0 missing_nums = missing_runs or all( [ - len([num for num in self.nums[run] if not num in saved_nums[run]]) + len([num for num in self.nums[run] if num not in saved_nums[run]]) > 0 for run in self.nums if run in saved_nums @@ -182,18 +193,18 @@ class Comparator(Aggregator, HDF5Container): def get_attr(self, attr_name, run, num, node_name="/", arg=None): pp = self.pp[run][num] - if not arg is None: + 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] - if not arg is None: + if arg is not None: name = name + "_" + str(arg) return pp.get_value(name) def get_global(self, node_name, run, num, arg=None, unload_cells=False): - if not arg is None: + if arg is not None: node_name = node_name + "_" + str(arg) pp = self.pp[run][num] if unload_cells: @@ -276,7 +287,6 @@ class Comparator(Aggregator, HDF5Container): return series def _from_log(self, keys, extractor): - nums = self.nums # Initialize series series = {} @@ -318,15 +328,15 @@ class Comparator(Aggregator, HDF5Container): for run in self.runs: # Surface of the box in pc^2 info = self.pp[run][self.nums[run][0]].info - surface = (info["unit_length"].express(cst.pc)) ** 2 + surface = (info["unit_length"].express(U.pc)) ** 2 # WARNING : We do not multiply by boxlen since already done in 'unit_length' (pymses) time = self.save.get_node("/series/sinks_from_log/time/" + run).read() - time = time * time_unit.express(cst.year) + time = time * time_unit.express(U.year) mass_sink = self.save.get_node( "/series/sinks_from_log/mass_sink/" + run ).read() - mass_sink = mass_sink * mass_unit.express(cst.Msun) + mass_sink = mass_sink * mass_unit.express(U.Msun) if avg_window is None: shift = 1 @@ -350,11 +360,11 @@ class Comparator(Aggregator, HDF5Container): for run in self.runs: # Surface of the box in pc^2 info = self.pp[run][self.nums[run][0]].info - surface = (info["unit_length"].express(cst.pc)) ** 2 + surface = (info["unit_length"].express(U.pc)) ** 2 mass_sink = self.save.get_node( "/series/sinks_from_log/mass_sink/" + run ).read() - mass_sink = mass_sink * mass_unit.express(cst.Msun) + mass_sink = mass_sink * mass_unit.express(U.Msun) ssm[run] = mass_sink / surface @@ -402,7 +412,7 @@ class Comparator(Aggregator, HDF5Container): glob_group="/globals", subarray_name=None, unload_cells=True, - unit=cst.none, + unit=U.none, description="", ): @@ -472,7 +482,7 @@ class Comparator(Aggregator, HDF5Container): ) def def_rules(self): - averageables = ["coldens", "rho", "T", "Q"] + self.rules = { # Read from log "sinks_from_log": Rule( @@ -483,7 +493,7 @@ class Comparator(Aggregator, HDF5Container): self._extract_sinks_from_log, ), group="/series", - unit={"time": "unit_time", "mass_sink": cst.Msun, "nb_sink": cst.none}, + unit={"time": "unit_time", "mass_sink": U.Msun, "nb_sink": U.none}, description={ "time": "Time", "mass_sink": "Total mass of stars", @@ -494,7 +504,7 @@ class Comparator(Aggregator, HDF5Container): self, self._ssfr_from_mass_sink, group="/series/sinks_from_log", - unit=cst.ssfr, + unit=U.ssfr, description="Instantaneous surfacic star formation rate", dependencies=["sinks_from_log"], ), @@ -502,7 +512,7 @@ class Comparator(Aggregator, HDF5Container): self, self._surfacic_sink_mass, group="/series/sinks_from_log", - unit=cst.Msun / cst.pc ** 2, + unit=U.Msun / U.pc ** 2, description="Surfacic sink mass", dependencies=["sinks_from_log"], ), @@ -510,7 +520,7 @@ class Comparator(Aggregator, HDF5Container): self, partial(self._from_log, ["time", "sfr"], self._extract_sfr_from_log), group="/series", - unit={"time": cst.year, "sfr": cst.ssfr}, + unit={"time": U.year, "sfr": U.ssfr}, description={ "time": "Time", "sfr": "Averaged surfacic star formation rate", @@ -527,7 +537,7 @@ class Comparator(Aggregator, HDF5Container): unit={ "time": "unit_time", "dt": "unit_time", - "turb_rms": cst.none, + "turb_rms": U.none, "turb_energy": { "unit_length": 3, "unit_velocity": 2, @@ -551,13 +561,13 @@ class Comparator(Aggregator, HDF5Container): group="/series", unit={ "time": "unit_time", - "step": cst.none, - "mcons": cst.none, - "econs": cst.none, - "epot": cst.none, # TODO find unit - "ekin": cst.none, - "eint": cst.none, - "emag": cst.none, + "step": U.none, + "mcons": U.none, + "econs": U.none, + "epot": U.none, # TODO find unit + "ekin": U.none, + "eint": U.none, + "emag": U.none, }, ), "turb_power": Rule( diff --git a/mypool.py b/mypool.py index ee26eaa..c6c2017 100644 --- a/mypool.py +++ b/mypool.py @@ -6,8 +6,6 @@ import multiprocessing # We must import this explicitly, it is not imported by the top-level # multiprocessing module. import multiprocessing.pool -import time -from random import randint class NoDaemonProcess(multiprocessing.Process): diff --git a/plotter.py b/plotter.py index 3576926..d4c5a10 100644 --- a/plotter.py +++ b/plotter.py @@ -9,9 +9,7 @@ This is the plotter module. """ import os -import sys from functools import partial - import matplotlib as mpl import numpy as np import tables @@ -27,10 +25,24 @@ if os.environ.get("DISPLAY", "") == "": mpl.use("Agg") import datetime -import pylab as P +import matplotlib.pyplot as plt import pspec_read -from comparator import * +from baseprocessor import Rule, BaseProcessor +from aggregator import Aggregator +from comparator import Comparator +from run_selector import RunSelector +from units import U, unit_str, convert_exp + + +from astrophysix.simdm.results import GenericResult +from astrophysix.simdm.experiment import ( + ParameterSetting, + ParameterVisibility, + Simulation, +) +from ramses_astrophysix import ramses + filetype_from_ext = {ext: ft for ft in FileType for ext in ft.extension_list} @@ -88,19 +100,18 @@ class Plotter(Aggregator, BaseProcessor): "beta": "$\\beta$", "beta_cool": "$\\beta$", "dens0": "$n_0$", - "coldens0": "$\Sigma_0$", + "coldens0": "$\\Sigma_0$", "sfr_avg_window": "window", "bx_bound": "$B_0$", - "levelmax": "$l_{\max}$", - "levelmin": "$l_{\min}$", + "levelmax": "$l_{\\max}$", + "levelmin": "$l_{\\min}$", "comp_frac": "$1 - \\zeta$", } # Conversion table from namelist values (from amses config file) into LaTex strings value_convert = { "sfr_avg_window": lambda x: "${:g}$ Myr".format(80 * x), - #'comp_frac' : lambda x: "${:g}$".format(1 - x), - "bx_bound": lambda x: "${:g}$ $\mu G$".format(5.267501272979475 * x), + "bx_bound": lambda x: "${:g}$ $\\mu G$".format(5.267501272979475 * x), } def __init__( @@ -112,7 +123,7 @@ class Plotter(Aggregator, BaseProcessor): pp_params=None, selector=None, tag=None, - unit_time=cst.year, + unit_time=U.year, **kwargs, ): @@ -124,10 +135,12 @@ class Plotter(Aggregator, BaseProcessor): path : path to the main folder of the simulations (ex '~/simus/myproject') in_runs : list of the runs to consider (ex ['run1', 'run2']) - in_nums : list or dict of the outputs numbers to consider (ex [3, 5] or {'run1' : [3, 5], 'run2' : [4, 6]) + 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. - selector : Existing instance of RunSelector, that selects runs and outputs. If set, in_runs and in_nums will be ignored + 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. """ @@ -216,7 +229,7 @@ class Plotter(Aggregator, BaseProcessor): dep, dep_arg, overwrite, overwrite_dep=self.overwrite_dep ) if result is not None: - self.just_done.append(done) + self.just_done.append(result) else: super(Plotter, self)._not_self_dep(name, dep, dep_arg, overwrite, **kwargs) @@ -238,7 +251,7 @@ class Plotter(Aggregator, BaseProcessor): """ # Set full name according to argument - if not arg is None: + if arg is not None: name_full = ( name + "_" @@ -286,7 +299,7 @@ class Plotter(Aggregator, BaseProcessor): real_ax = ax[i] except TypeError as e: if ax is None: - fig, real_ax = P.subplots(1, 1) + fig, real_ax = plt.subplots(1, 1) elif not_array_error(e): real_ax = ax else: @@ -341,15 +354,15 @@ class Plotter(Aggregator, BaseProcessor): """ Once all dependencies are met, actually process the rule """ - P.sca(ax) + plt.sca(ax) if self._needs_computation(overwrite, plot_filename): plot_info = rule.plot(save, arg, **kwargs) if not self.pp_params.out.interactive: - P.tight_layout(pad=1) + plt.tight_layout(pad=1) if self.pp_params.out.save: - P.savefig(plot_filename) + plt.savefig(plot_filename) self._log("{} plotted".format(plot_filename), "SUCCESS") else: self._log( @@ -357,7 +370,7 @@ class Plotter(Aggregator, BaseProcessor): ) if not self.pp_params.out.interactive: - P.close() + plt.close() return plot_info else: self._log("Plot {} is already done, skipping...".format(plot_filename)) @@ -372,9 +385,9 @@ class Plotter(Aggregator, BaseProcessor): if not self.pp_params.out.tag == "": tag_name = "_" + tag_name - if not run is None and not num is None: + if run is not None and num is not None: fmt = "{out}/{run}/{name}{tag}_{run}_{num:05}{ext}" - elif not run is None: + elif run is not None: fmt = "{out}/{run}/{name}{tag}_{run}{ext}" else: fmt = "{out}/{name}{tag}{ext}" @@ -382,7 +395,7 @@ class Plotter(Aggregator, BaseProcessor): fmt = self.pp_params.out.fmt nml = None - if not run is None: + if run is not None: nml = self.comp.namelist[run] return fmt.format( @@ -423,12 +436,12 @@ class Plotter(Aggregator, BaseProcessor): label_run = r"{}".format(self.save.root._v_attrs.attrs[run].label) else: label_run = run - elif not nml_key is None: + elif nml_key is not None: if not type(nml_key) == list: nml_key = [nml_key] label_run = ", ".join(map(get_label_nml, nml_key)) - if not label is None: + if label is not None: label_run = label + " (" + label_run + ")" else: label_run = label @@ -451,7 +464,7 @@ class Plotter(Aggregator, BaseProcessor): if "unit" in node._v_attrs: unit_old = node._v_attrs.unit else: - unit_old = cst.none + unit_old = U.none if unit is None: unit = unit_old @@ -466,7 +479,7 @@ class Plotter(Aggregator, BaseProcessor): return label, unit_old, unit - def _snapshot_title(self, run, node, title, nml_key, put_time, unit_time=cst.Myr): + def _snapshot_title(self, run, node, title, nml_key, put_time, unit_time=U.Myr): title = self._label_run(run, node, title, nml_key) if put_time: @@ -497,9 +510,9 @@ class Plotter(Aggregator, BaseProcessor): put_title=True, nml_key=None, put_time=True, - unit_time=cst.Myr, + unit_time=U.Myr, put_units=True, - unit_space=cst.pc, + unit_space=U.pc, cmap="plasma", norm="log", put_cbar=True, @@ -535,14 +548,14 @@ class Plotter(Aggregator, BaseProcessor): elif norm == "linear": norm = mpl.colors.NoNorm() - if autoscale and not norm is None: + if autoscale and norm is not None: norm.autoscale(dmap) - im = P.imshow( + im = plt.imshow( dmap, extent=im_extent, origin="lower", norm=norm, cmap=cmap, **kwargs ) - P.locator_params(axis="both", nbins=self.pp_params.plot.ntick) + plt.locator_params(axis="both", nbins=self.pp_params.plot.ntick) if xlabel is None: xlabel = self._ax_title[ax_h] @@ -551,19 +564,19 @@ class Plotter(Aggregator, BaseProcessor): if put_units: xlabel = xlabel + unit_str(unit_space) ylabel = ylabel + unit_str(unit_space) - P.xlabel(xlabel) - P.ylabel(ylabel) + plt.xlabel(xlabel) + plt.ylabel(ylabel) try: - cbar = P.colorbar(im, cax=P.gca().cax) + cbar = plt.colorbar(im, cax=plt.gca().cax) except AttributeError: - cbar = P.colorbar() + cbar = plt.colorbar() if put_title: title = self._snapshot_title(run, node, title, nml_key, put_time, unit_time) - P.title(title) + plt.title(title) - if not label is None: + if label is not None: cbar.set_label(label) for i, plot_overlay in enumerate(overlays): @@ -627,7 +640,7 @@ class Plotter(Aggregator, BaseProcessor): ) lw[lvl_array < lvl_th] = 1.0 - cont = P.contour( + cont = plt.contour( map_contour, extent=im_extent, origin="lower", @@ -639,7 +652,7 @@ class Plotter(Aggregator, BaseProcessor): lvls = np.array(cont.levels) + lvl_offset cont.levels = lvls - P.clabel( + plt.clabel( cont, lvls[np.array(lvls) < lvl_max_lbl], inline=1, @@ -663,13 +676,11 @@ class Plotter(Aggregator, BaseProcessor): ) def _overlay_speed( - self, ax_los, im_extent, unit=cst.km_s, unit_coeff=1.0, key_v=None, **kwargs + self, ax_los, im_extent, unit=U.km_s, unit_coeff=1.0, key_v=None, **kwargs ): """ Add an overlay : velocity vector field """ - ax_h = self._axes_h[ax_los] - ax_v = self._axes_v[ax_los] dmap_vh_node = self.save.get_node("/maps/speed_h_{}".format(ax_los)) dmap_vh = dmap_vh_node.read() dmap_vv = self.save.get_node("/maps/speed_v_{}".format(ax_los)).read() @@ -701,11 +712,11 @@ class Plotter(Aggregator, BaseProcessor): max_v = np.max(norm_v) min_v = np.min(norm_v) - Q = P.quiver(hh, vv, map_vh_red, map_vv_red, units="width", **kwargs) + Q = plt.quiver(hh, vv, map_vh_red, map_vv_red, units="width", **kwargs) if key_v is None: key_v = (max_v + min_v) / 2.0 - P.quiverkey( + plt.quiverkey( Q, 0.6, 0.98, @@ -719,8 +730,6 @@ class Plotter(Aggregator, BaseProcessor): """ Add an overlay : magnetic streamlines """ - ax_h = self._axes_h[ax_los] - ax_v = self._axes_v[ax_los] dmap_Bh_node = self.save.get_node("/maps/B_h_{}".format(ax_los)) dmap_Bh = dmap_Bh_node.read() dmap_Bv = self.save.get_node("/maps/B_v_{}".format(ax_los)).read() @@ -744,7 +753,7 @@ class Plotter(Aggregator, BaseProcessor): ) * lbox hh, vv = np.meshgrid(vec_h, vec_v) - P.streamplot(hh, vv, map_Bh_red, map_Bv_red, **kwargs) + plt.streamplot(hh, vv, map_Bh_red, map_Bv_red, **kwargs) def _plot_radial( self, @@ -760,7 +769,7 @@ class Plotter(Aggregator, BaseProcessor): nml_key=None, put_title=True, put_time=True, - unit_time=cst.Myr, + unit_time=U.Myr, **kwargs, ): """ @@ -775,23 +784,23 @@ class Plotter(Aggregator, BaseProcessor): if ytransform is not None: mean_bin = ytransform(mean_bin) - P.xlabel(r"$r$") + plt.xlabel(r"$r$") if xlog: - P.xscale("log") + plt.xscale("log") if ylog: - P.yscale("log") + plt.yscale("log") - if not ylabel is None: - P.ylabel(ylabel) + if ylabel is not None: + plt.ylabel(ylabel) title = self._snapshot_title(run, node, title, nml_key, put_time, unit_time) if put_title: - P.title(title) - if label == None: + plt.title(title) + if label is None: label = title - P.plot(bin_centers, mean_bin, label=label, **kwargs) + plt.plot(bin_centers, mean_bin, label=label, **kwargs) def _plot_hist( self, @@ -808,11 +817,11 @@ class Plotter(Aggregator, BaseProcessor): title=None, nml_key=None, put_time=True, - unit_time=cst.Myr, + unit_time=U.Myr, xlog=None, ylog=False, kind="bar", - ylabel="$\mathcal{P}$", + ylabel="$\\mathcal{P}$", color=None, colors=None, nml_color=None, @@ -824,7 +833,7 @@ class Plotter(Aggregator, BaseProcessor): Plot an histogram (PDF, etc ...) """ # Get node - if not ax_los is None: + if ax_los is not None: name = name + "_" + ax_los node = self.save.get_node(group + name) if xlog is None: @@ -853,12 +862,12 @@ class Plotter(Aggregator, BaseProcessor): # Set title title = self._snapshot_title(run, node, title, nml_key, put_time, unit_time) if put_title: - P.title(title) - if label == None: + plt.title(title) + if label is None: label = title # Set colors - if color is None and not colors is None: + if color is None and colors is not None: if nml_color is None: color = colors[run] else: @@ -870,25 +879,27 @@ class Plotter(Aggregator, BaseProcessor): # Actual plot if kind == "bar": - P.bar(centers, values, width, log=ylog, color=color, label=label, **kwargs) + plt.bar( + centers, values, width, log=ylog, color=color, label=label, **kwargs + ) elif kind == "step": if ylog: - P.yscale("log") - P.step(centers, values, where="mid", color=color, label=label, **kwargs) + plt.yscale("log") + plt.step(centers, values, where="mid", color=color, label=label, **kwargs) else: raise ValueError("kind must be 'bar' or 'step'") # put labels - if not label is None: - P.xlabel(xlabel) - if not ylabel is None: - P.ylabel(ylabel) + if label is not None: + plt.xlabel(xlabel) + if ylabel is not None: + plt.ylabel(ylabel) # Also diplay fit, previously saved if ax_los is not None and "/hist/fit_" + name + "_" + ax_los in self.save: slope = node.attrs.slope origin = node.attrs.origin - P.plot( + plt.plot( centers, 10 ** (slope * centers + origin), "--", @@ -896,7 +907,7 @@ class Plotter(Aggregator, BaseProcessor): color="orange", ) # or a new one - if not fit is None: + if fit is not None: self._overlay_fit( centers, values, kind=fit, ls="--", lw=1.5, label=fitlabel ) @@ -912,7 +923,7 @@ class Plotter(Aggregator, BaseProcessor): xaxis_label=xlabel, yaxis_label=ylabel, xaxis_unit=unit, - yaxis_unit=cst.none, + yaxis_unit=U.none, plot_title=title, ) @@ -938,7 +949,7 @@ class Plotter(Aggregator, BaseProcessor): sigma_err=2.0, grid=False, put_time=False, - unit_time=cst.Myr, + unit_time=U.Myr, colors=None, nml_color=None, legend=None, @@ -951,7 +962,7 @@ class Plotter(Aggregator, BaseProcessor): """ # Get proper hdf5 names - if not node_arg is None: + if node_arg is not None: name_x, name_y = name_x + "_" + node_arg, name_y + "_" + node_arg # Get hdf5 nodes @@ -1031,15 +1042,15 @@ class Plotter(Aggregator, BaseProcessor): if smooth > 0: y = gaussian_filter1d(y, sigma=smooth) - if not run is None: + if run is not None: label = self._label_run(run, node_y, label, nml_key) # Look if special colors method is used if colors is None: if yerr is None: - (base_line,) = P.plot(x, y, label=label, **kwargs) + (base_line,) = plt.plot(x, y, label=label, **kwargs) else: - base_line, _, _ = P.errorbar(x, y, yerr=yerr, label=label, **kwargs) + base_line, _, _ = plt.errorbar(x, y, yerr=yerr, label=label, **kwargs) else: if nml_color is None: color = colors[run] @@ -1055,21 +1066,21 @@ class Plotter(Aggregator, BaseProcessor): except: color = colors(nml) if yerr is None: - (base_line,) = P.plot(x, y, label=label, color=color, **kwargs) + (base_line,) = plt.plot(x, y, label=label, color=color, **kwargs) else: - base_line, _, _ = P.errorbar( + base_line, _, _ = plt.errorbar( x, y, yerr=yerr, color=color, label=label, **kwargs ) # Ax decorations - P.xlabel(xlabel) - P.ylabel(ylabel) + plt.xlabel(xlabel) + plt.ylabel(ylabel) if grid: - P.grid() + plt.grid() if legend: - P.legend() + plt.legend() - if not fit is None: + if fit is not None: self._overlay_fit( x, y, @@ -1121,7 +1132,7 @@ class Plotter(Aggregator, BaseProcessor): ) if label is None: label = r"Linear fit with slope ${:.3g}$".format(a) - P.plot(x, a * x + b, label=label, **kwargs) + plt.plot(x, a * x + b, label=label, **kwargs) elif kind == "power_law": if yerr is None: (a, b, rho, _map_rule, stderr) = linregress(np.log10(x), np.log10(y)) @@ -1131,8 +1142,13 @@ class Plotter(Aggregator, BaseProcessor): ) ) else: - fitfunc = lambda p, x: p[0] + p[1] * x - errfunc = lambda p, x, y, err: (y - fitfunc(p, x)) / err + + def fitfunc(p, x): + return p[0] + p[1] * x + + def errfunc(p, x, y, err): + return (y - fitfunc(p, x)) / err + pinit = [1.0, -1.0] out = optimize.leastsq( errfunc, @@ -1151,27 +1167,27 @@ class Plotter(Aggregator, BaseProcessor): ) if label is None: label = r"Power-law fit with index {:.1f}".format(a) - P.plot(x, (10 ** b) * x ** a, label=label, **kwargs) + plt.plot(x, (10 ** b) * x ** a, label=label, **kwargs) def overlay_kennicutt(self, n0, step): """ Add an overlay : Kennicutt mass accretion """ - P.grid(False) - ylim = P.ylim() - (tmin, tmax) = P.xlim() + plt.grid(False) + ylim = plt.ylim() + (tmin, tmax) = plt.xlim() tmax = tmax + 20 - ymax = P.ylim()[1] + ymax = plt.ylim()[1] ssfr_sun = 2.5e-9 ssfr_ken = ssfr_sun * n0 ** 1.4 - coeff = ssfr_ken * 1e6 * (self.comp.info["unit_length"].express(cst.pc)) ** 2 + coeff = ssfr_ken * 1e6 * (self.comp.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) - P.plot(t + i, t * coeff, ls="--", lw=0.9, color="grey") - P.plot(t + tmin, (t + i - tmin) * coeff, ls="--", lw=0.9, color="grey") - P.xlim(tmin, tmax) - P.ylim(ylim) + plt.plot(t + i, t * coeff, ls="--", lw=0.9, color="grey") + plt.plot(t + tmin, (t + i - tmin) * coeff, ls="--", lw=0.9, color="grey") + plt.xlim(tmin, tmax) + plt.ylim(ylim) def _gen_from_log(self, logrule, name, description="Generated"): self.rules[name] = PlotRule( @@ -1180,7 +1196,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot, "/series/" + logrule + "/time", "/series/" + logrule + "/" + name, - xunit=cst.Myr, + xunit=U.Myr, ), description=description, kind="series", @@ -1202,7 +1218,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot_map, "coldens", label=r"$\Sigma$", - # unit=cst.coldens + # unit=U.coldens ), "Column density map", dependencies=["coldens"], @@ -1235,7 +1251,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot_map, "vphi", label=r"$v_\phi$", - # unit=cst.km_s + # unit=U.km_s ), "Azimuthal speed", dependencies=["vphi"], @@ -1246,7 +1262,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot_map, "vr", label=r"$v_r$", - # unit=cst.km_s + # unit=U.km_s ), "Radial speed", dependencies=["vr"], @@ -1257,7 +1273,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot_map, "rho", label=r"$\rho$", - # unit=cst.Msun_pc3 + # unit=U.Msun_pc3 ), "Density slice at s = 0, with s = x, y or z.", dependencies=["rho"], @@ -1268,7 +1284,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot_map, "coldens", label=r"$\Sigma$", - unit=cst.coldens, + unit=U.coldens, overlays=[self._overlay_levels], ), "Column density with level overlay", @@ -1280,7 +1296,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot_map, "rho", label=r"$\rho$", - unit=cst.Msun_pc3, + unit=U.Msun_pc3, overlays=[self._overlay_speed], ), "Density slice with speed overlay", @@ -1292,7 +1308,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot_map, "rho", label=r"$\rho$", - unit=cst.Msun_pc3, + unit=U.Msun_pc3, overlays=[self._overlay_B], ), "Density slice with magnetic field overlay", @@ -1304,7 +1320,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot_map, "rho", label=r"$\rho$", - unit=cst.Msun_pc3, + unit=U.Msun_pc3, overlays=[self._overlay_B, self._overlay_speed], ), "Density slice with magnetic field and velocity overlay", @@ -1382,7 +1398,7 @@ class Plotter(Aggregator, BaseProcessor): "B_int": PlotRule( self, partial( - self._plot_map, "B_int", label=r"$\mid \mathrm{B} \mid$", unit=cst.T + self._plot_map, "B_int", label=r"$\mid \mathrm{B} \mid$", unit=U.T ), "Magnetic intensity map", dependencies=["B_int"], @@ -1450,8 +1466,8 @@ class Plotter(Aggregator, BaseProcessor): self._plot, "/series/sinks_from_log/time", "/series/sinks_from_log/mass_sink", - xunit=cst.Myr, - yunit=cst.Msun, + xunit=U.Myr, + yunit=U.Msun, ), "Mass of the sinks as a function of time", kind="series", @@ -1463,8 +1479,8 @@ class Plotter(Aggregator, BaseProcessor): self._plot, "/series/sinks_from_log/time", "/series/sinks_from_log/ssm", - xunit=cst.Myr, - yunit=cst.Msun / cst.pc ** 2, + xunit=U.Myr, + yunit=U.Msun / U.pc ** 2, ), "Mass of the sinks as a function of time divided by surface", kind="series", @@ -1477,8 +1493,8 @@ class Plotter(Aggregator, BaseProcessor): "/series/sfr_from_log/time", "/series/sfr_from_log/sfr", ylabel="Averaged surfacic SFR", - xunit=cst.Myr, - yunit=cst.ssfr, + xunit=U.Myr, + yunit=U.ssfr, ), kind="series", dependencies=["sfr_from_log"], @@ -1490,8 +1506,8 @@ class Plotter(Aggregator, BaseProcessor): "/series/sinks_from_log/time", "/series/sinks_from_log/issfr", ylabel="Surfacic SFR", - xunit=cst.Myr, - yunit=cst.ssfr, + xunit=U.Myr, + yunit=U.ssfr, ), kind="series", dependencies=["issfr"], @@ -1502,7 +1518,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot, "/series/rms_from_log/time", "/series/rms_from_log/turb_rms", - xunit=cst.Myr, + xunit=U.Myr, ), "Turbulent RMS", kind="series", @@ -1514,7 +1530,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot, "/series/rms_from_log/time", "/series/rms_from_log/turb_energy", - xunit=cst.Myr, + xunit=U.Myr, ), "Turbulent energy", kind="series", @@ -1526,7 +1542,7 @@ class Plotter(Aggregator, BaseProcessor): self._plot, "/series/rms_from_log/time", "/series/rms_from_log/turb_power", - xunit=cst.Myr, + xunit=U.Myr, ), "Turbulent power", kind="series", @@ -1539,8 +1555,8 @@ class Plotter(Aggregator, BaseProcessor): "/series/time", "/series/time_sigma", ylabel="$\\sigma$", - xunit=cst.Myr, - yunit=cst.km_s, + xunit=U.Myr, + yunit=U.km_s, ), "Velocity dispersion", kind="series", @@ -1552,8 +1568,8 @@ class Plotter(Aggregator, BaseProcessor): self._plot, "/series/time", "/series/time_mwa_B_int", - xunit=cst.Myr, - yunit=cst.uG, + xunit=U.Myr, + yunit=U.uG, ), "Magnetic intensity average", kind="series", @@ -1565,8 +1581,8 @@ class Plotter(Aggregator, BaseProcessor): self._plot, "/series/time", "/series/time_mass", - xunit=cst.Myr, - yunit=cst.Msun, + xunit=U.Myr, + yunit=U.Msun, ), "Total mass in the box", kind="series", @@ -1578,8 +1594,8 @@ class Plotter(Aggregator, BaseProcessor): self._plot, "/series/time", "/series/time_max_fluct_coldens_z", - ylabel="$\\max(\Sigma/\overline{\Sigma})$", - xunit=cst.Myr, + ylabel="$\\max(\\Sigma/\\overline{\\Sigma})$", + xunit=U.Myr, ), "Maximal fluctuation of the column density against time", kind="series", diff --git a/postprocessor.py b/postprocessor.py index 3632e3d..135c642 100644 --- a/postprocessor.py +++ b/postprocessor.py @@ -1,15 +1,39 @@ # coding: utf-8 -import pickle +import numpy as np +import tables +import pickle import astropy.units as u import pandas as pd -import pymses.utils.regions as reg -from fil_finder import FilFinder2D -from pymses.filters import RegionFilter + from skimage.morphology import medial_axis +import os +from functools import partial +from scipy.stats import linregress + +from astrophysix.simdm.results import Snapshot + +import pymses +import pymses.utils.regions as reg +from pymses.analysis import ( + Camera, + FractionOperator, + MaxLevelOperator, + ScalarOperator, + raytracing, + slicing, + splatting, +) +from pymses.filters import CellsToPoints, RegionFilter + +from fil_finder import FilFinder2D import pspec_new -from baseprocessor import * + + +from units import U +from baseprocessor import HDF5Container, Rule, norm_getter, simple_getter, vect_getter + # Getters @@ -167,7 +191,7 @@ class PostProcessor(HDF5Container): path_out=None, pp_params=None, tag=None, - unit_time=cst.year, + unit_time=U.year, ): """ Creates the basic structures needed for the outputs @@ -246,21 +270,25 @@ class PostProcessor(HDF5Container): if self.pp_params.pymses.filter: center = (self.max_coords + self.min_coords) / 2.0 - im_extent = [ - self.min_coords[0], - self.max_coords[0], - self.min_coords[1], - self.max_coords[1], - ] + 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 = [ - (-self._radius + center[0]), - (self._radius + center[0]), - (-self._radius + center[1]), - (self._radius + center[1]), - ] + 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 @@ -274,7 +302,7 @@ class PostProcessor(HDF5Container): self.save.root._v_attrs.unit_length = self.info["unit_length"] self.save.root._v_attrs.time = time - if not "/maps" in self.save: + 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 @@ -283,7 +311,6 @@ class PostProcessor(HDF5Container): # Initialize cameras self._cam = {} for ax_los in self._ax_nb: # los = line of sight - ax_h = self._axes_h[ax_los] ax_v = self._axes_v[ax_los] self._cam[ax_los] = Camera( @@ -306,13 +333,13 @@ class PostProcessor(HDF5Container): else: self.fil = None + time_in_right_unit = self.info["time"] * self.info["unit_time"].express( + unit_time + ) self.snapshot = Snapshot( name=str(self.num), description="", - time=( - self.info["time"] * self.info["unit_time"].express(unit_time), - unit_time, - ), + time=(time_in_right_unit, unit_time), directory_path=self.path, data_reference="OUTPUT_{}".format(self.num), ) @@ -323,7 +350,7 @@ class PostProcessor(HDF5Container): """ Load all cells from the source file in the memory. Cells will be accessible trough self.cells - (/!\ Long and memory heavy) + (Long and memory heavy) """ if not self.cells_loaded: if os.path.exists(self.cells_filename): @@ -424,7 +451,7 @@ class PostProcessor(HDF5Container): """ Azimuthal velocity """ return self.oct_getter_vect_phi(dset, "vel") - def _slice(self, getter, ax_los="z", z=0.0, unit=cst.none): + def _slice(self, getter, ax_los="z", z=0.0, unit=U.none): """ Slice process function. Return a slice of the source box. @@ -440,7 +467,7 @@ class PostProcessor(HDF5Container): z : float Coordinate of the slice on the ax_los axis - unit : cst.Unit + unit : U.Unit Unit of the resulting dataset Returns @@ -452,9 +479,7 @@ class PostProcessor(HDF5Container): datamap = slicing.SliceMap(self._amr, self._cam[ax_los], op, z=z) return datamap.map.T - def _ax_avg( - self, getter, ax_los, unit=cst.none, mass_weighted=True, surf_qty=False - ): + def _ax_avg(self, getter, ax_los, unit=U.none, mass_weighted=True, surf_qty=False): """ Map of the average of a quantity (given by getter) along an axis (ax_los) Returns 2D array if getter returns a scalar quantity @@ -505,7 +530,7 @@ class PostProcessor(HDF5Container): self.load_cells() return np.sort(np.unique(self.cells["pos"][:, axis])) - def _plane_avg_uniform(self, getter, axis, unit=cst.none, mass_weighted=False): + def _plane_avg_uniform(self, getter, axis, unit=U.none, mass_weighted=False): """ Profile of the average of a quantity (given by getter) perpendicular to an axis WARNING : This version only works on an uniform grid, need of a box version for AMR @@ -622,17 +647,17 @@ class PostProcessor(HDF5Container): """ self.load_cells() mean_speed = self.save.get_node("/globals/mwa_speed").read() - mean_speed = mean_speed * self.info["unit_velocity"].express(cst.km_s) + mean_speed = mean_speed * self.info["unit_velocity"].express(U.km_s) vel_fluct = (self.cells)["vel"] * self.info["unit_velocity"].express( - cst.km_s + U.km_s ) - mean_speed B_norm = getter_B_int(self.cells) - B_norm = B_norm * self.info["unit_mag"].express(cst.T) + B_norm = B_norm * self.info["unit_mag"].express(U.T) v_norm = np.sqrt( np.sum((vel_fluct * 10 ** (3)) ** 2, axis=1) ) # v_norm [m/s] et vel_fluct [km/s] rho = getter_rho(self.cells) - rho_kg_m3 = rho * self.info["unit_density"].express(cst.kg_m3) + rho_kg_m3 = rho * self.info["unit_density"].express(U.kg_m3) eb = 0.5 * (B_norm) ** 2 / (4 * np.pi * 10 ** (-7)) # mettre le bon mu ek = 0.5 * v_norm ** 2 * rho_kg_m3 rapport = ek / eb @@ -1024,8 +1049,7 @@ class PostProcessor(HDF5Container): return alpha alpha_f = ( - self._ax_avg(getter_alpha_num, "z", unit=cst.none, mass_weighted=True) - / T_avg + self._ax_avg(getter_alpha_num, "z", unit=U.none, mass_weighted=True) / T_avg ) # alpha @@ -1047,7 +1071,7 @@ class PostProcessor(HDF5Container): gphi = self.oct_getter_vect_phi(dset, "g") return gr * gphi / (4 * np.pi * self.G) - alpha_g = self._ax_avg(getter_alpha_grav, "z", unit=cst.none, surf_qty=True) / ( + alpha_g = self._ax_avg(getter_alpha_grav, "z", unit=U.none, surf_qty=True) / ( coldens * T_avg ) @@ -1166,10 +1190,6 @@ class PostProcessor(HDF5Container): GM = self.G * self.pp_params.disk.mass_star # Mass parameter - # Get mask for filaments - fil = self.fil - mask_fil = np.asarray(fil.mask.copy(), dtype=bool) - # Find center of filaments i_center, j_center = self._filaments_center() @@ -1283,7 +1303,7 @@ class PostProcessor(HDF5Container): self._alpha_disk, "Map of the Shakura&Sunaev alpha parameter for disks", "/maps", - unit=cst.none, + unit=U.none, dependencies=[ "avg_map_rho_avg", "avg_map_T_mwavg", @@ -1297,7 +1317,7 @@ class PostProcessor(HDF5Container): "Map of the graviational contrib to\ Shakura&Sunaev alpha parameter for disks", "/maps", - unit=cst.none, + unit=U.none, dependencies=["avg_map_coldens", "avg_map_T_mwavg"], ), "rho": Rule( @@ -1373,9 +1393,9 @@ class PostProcessor(HDF5Container): self._sinks, group="/datasets", unit={ - "Id": cst.none, - "M": cst.Msun, - "dmf": cst.Msun, + "Id": U.none, + "M": U.Msun, + "dmf": U.Msun, "x": "", "y": "", "z": "", @@ -1388,9 +1408,9 @@ class PostProcessor(HDF5Container): "lz": "|l|", "acc_rate": "[Msol/y]", "acc_lum": "[Lsol]", - "age": cst.year, + "age": U.year, "int_lum": "[Lsol]", - "Teff": cst.K, + "Teff": U.K, }, ), "pspec": Rule(self, self._pspec, "Power spectrum", "/hdf5"), @@ -1473,7 +1493,7 @@ class PostProcessor(HDF5Container): "Global cos fluctuation-PDF", "/hist", dependencies=["mwa_speed"], - unit=cst.none, + unit=U.none, ), "Brho": Rule( self, @@ -1488,7 +1508,7 @@ class PostProcessor(HDF5Container): "Average of Ek/Eb as a function of rho", "/datasets", dependencies=["mwa_speed"], - unit={"rho": self.info["unit_density"], "Ek_Eb_rho": cst.none}, + unit={"rho": self.info["unit_density"], "Ek_Eb_rho": U.none}, ), # Profiles "axis": Rule( diff --git a/pp_params.py b/pp_params.py index 369186b..6be735d 100644 --- a/pp_params.py +++ b/pp_params.py @@ -4,7 +4,6 @@ import os import re import munch -import numpy as np import yaml _dir_path = os.path.dirname(os.path.realpath(__file__)) diff --git a/run_selector.py b/run_selector.py index 644065a..702795a 100644 --- a/run_selector.py +++ b/run_selector.py @@ -5,10 +5,12 @@ import glob import os from functools import partial +import numpy as np +import yaml import f90nml -from pp_params import * +from pp_params import default_params class NamelistRecursive: @@ -247,7 +249,7 @@ class RunSelector: runs = self.nml_select(runs, filter_nml) # Sort by the value in the namelist of sort_run_by - if not sort_run_by is None: + if sort_run_by is not None: if type(sort_run_by) == str: sort_run_by = [sort_run_by] for nml_key in reversed(sort_run_by): diff --git a/units.py b/units.py index cafe8a8..8a2fed3 100644 --- a/units.py +++ b/units.py @@ -1,13 +1,13 @@ # coding: utf-8 -import astrophysix.units as cst +import astrophysix.units as U -create_unit = cst.Unit.create_unit +create_unit = U.Unit.create_unit def parse_exp_unit(u): splitted = u.split("^") - name_u = cst.Unit.from_name(splitted[0]).latex.replace("text", "math") + name_u = U.Unit.from_name(splitted[0]).latex.replace("text", "math") exp = "" if len(splitted) > 1: exp = "^{" + str(splitted[1]) + "}" @@ -39,7 +39,7 @@ def unit_str(unit, base=None, prefix="", format=" [{unit}]"): prefix : str to put befor the unit format : str with the {unit} key, to put external decoration """ - if unit == cst.none: + if unit == U.none: return "" elif not base is None: coeff = unit.express(base) @@ -64,33 +64,31 @@ def unit_str(unit, base=None, prefix="", format=" [{unit}]"): return format.format(unit=u_str) -cst.coldens = create_unit( - "Msun.pc^-2", base_unit=cst.Msun / cst.pc ** 2, descr="Column density" +U.coldens = create_unit( + "Msun.pc^-2", base_unit=U.Msun / U.pc ** 2, descr="Column density" ) -cst.km_s = create_unit("km.s^-1", base_unit=cst.km / cst.s, descr="Speed") +U.km_s = create_unit("km.s^-1", base_unit=U.km / U.s, descr="Speed") -cst.Msun_pc3 = create_unit( - "Msun.pc^-3", base_unit=cst.Msun / cst.pc ** 3, descr="Density" -) +U.Msun_pc3 = create_unit("Msun.pc^-3", base_unit=U.Msun / U.pc ** 3, descr="Density") -cst.kg_m3 = create_unit("kg.m^-3", base_unit=cst.kg / cst.m ** 3, descr="Density") +U.kg_m3 = create_unit("kg.m^-3", base_unit=U.kg / U.m ** 3, descr="Density") -cst.ssfr = create_unit( +U.ssfr = create_unit( "Msun.year^-1.pc^-2", - base_unit=cst.Msun / cst.year / cst.pc ** 2, + base_unit=U.Msun / U.year / U.pc ** 2, descr="Surfacic SFR", ) # latex='M$_{\odot}$.yr$^{-1}$.pc$^{-2}$') -cst.ssfrG = create_unit( +U.ssfrG = create_unit( "Msun.Gyr^-1.pc^-2", - base_unit=1e-9 * cst.Msun / cst.year / cst.pc ** 2, + base_unit=1e-9 * U.Msun / U.year / U.pc ** 2, descr="Surfacic SFR", latex="\mathrm{M}_{\odot}.\mathrm{Gyr}^{-1}.\mathrm{pc}^{-2}", ) -cst.uG = create_unit( - "μG", base_unit=1e-10 * cst.T, descr="Micro Gauss", latex="\\mu\\mathrm{G}" +U.uG = create_unit( + "μG", base_unit=1e-10 * U.T, descr="Micro Gauss", latex="\\mu\\mathrm{G}" )