diff --git a/disk_postprocess.py b/disk_postprocess.py index 14c6f86..78e7b70 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -5,11 +5,10 @@ import os import pymses import numpy as np import matplotlib as mpl - -if os.environ.get("DISPLAY", "") == "": - print("No display found. Using non-interactive Agg backend") - mpl.use("Agg") - +if os.environ.get('DISPLAY','') == '': + print('No display found. Using non-interactive Agg backend') + mpl.use('Agg') +from matplotlib.patches import Polygon import pylab as P import glob as glob @@ -27,6 +26,9 @@ from pymses.analysis import Camera, raytracing, slicing, splatting from pymses.filters import CellsToPoints from pymses.analysis import ScalarOperator, FractionOperator, MaxLevelOperator from scipy.stats import gaussian_kde +import matplotlib.patches as mpatches +from matplotlib.collections import PatchCollection + import module_extract as me @@ -119,22 +121,12 @@ def compute_image_data(path, num, radius=0.5, # Checking for existing files name_save = directory + '/maps_disk' + '_' + tag + '_' + format(num,'05') + '.save' if (len(glob.glob(name_save))==1 and not force): - return + return # Prepare saving data - if "T" in images and not "rho" in images: - images.append("rho") - if "Q" in images and not "coldens" in images: - images.append("coldens") - maps_disk = { - "time": time, - "im_extent": im_extent, - "center": center, - "radius": radius, - "lbox": lbox, - "images": images, - "axes_los": axes_los, - } + maps_disk = {'time' : time, 'im_extent' : im_extent, + 'center' : center, 'radius' : radius, 'lbox' : lbox, + 'images': images, 'axes_los': axes_los} # Prepare raytracing rho_op = ScalarOperator(lambda dset: dset["rho"], ro.info["unit_density"]) @@ -344,6 +336,11 @@ def plot_maps(path, num, ): continue + if (image == 'jeans'): + maps_disk['jeans_' + ax_los] = np.sqrt(np.pi * maps_disk['T_' + ax_los] / maps_disk['rho_' + ax_los]) + if (image == 'jeans_ratio'): + maps_disk['jeans_ratio_' + ax_los] = maps_disk['jeans_' + ax_los] * 2**(maps_disk['levels_' + ax_los]) + map_disk = maps_disk[image + '_' + ax_los] if image == 'Q' : @@ -353,6 +350,13 @@ def plot_maps(path, num, cmap='RdBu', norm=mpl.colors.LogNorm(), vmin=0.01, vmax=100.) + elif image == 'jeans_ratio' : + im = P.imshow(map_disk, + extent=im_extent, + origin='lower', + cmap='RdBu', + norm=mpl.colors.LogNorm(), + vmin=0.1, vmax=1000.) else: im = P.imshow( map_disk, @@ -396,7 +400,7 @@ def plot_maps(path, num, cont.levels = cont.levels + 1 P.clabel(cont, - levels_ar[levels_ar < 11][1:], + cont.levels[cont.levels < 11], inline=1, fontsize=8., fmt='%1d') elif image == 'rho': cbar.set_label(r'$\rho$ (code)') @@ -421,6 +425,8 @@ def plot_maps(path, num, cbar.set_label(r'$T (code)$') elif image == 'Q': cbar.set_label(r'$Q$') + elif image == 'jeans': + cbar.set_label(r'Jeans\'s lenght') else: cbar.set_label(image) @@ -788,7 +794,7 @@ def plot_disk_prop( P.plot(rad,abs(prop_disk['Q_kepl']),color='b',linewidth=2) # P.plot(rad, abs(prop_disk['Q_mean']) * np.ones(len(rad)), 'b:', linewidth=1) P.ylabel(r'$Q$') - title = title + title = title elif plot == 'H': P.plot(rad,abs(prop_disk['height']/ rad),color='b',linewidth=2) P.ylabel(r'H ratio') @@ -805,20 +811,67 @@ def plot_disk_prop( P.close() -def disk_pdf( - path, - num, - maps_disk, - pos_star=[1.0, 1.0], - force=False, - interactive=False, - nb_bin_hist=50, - tag="", - rad_min=0.075, - rad_max=0.3, - put_title=True, - do_speed=True, -): +def sigma_pdf(fluct_maps, mask_flat, put_title=False, title='', nb_bin_hist=50, tag=''): + # Sigma-PDF + dcoldens = np.log10(fluct_maps['coldens']).flatten() + + nb_cells = np.sum(mask_flat) + P.grid() + P.yscale('log') + P.ylim([0.5 / nb_cells, 1.]) + P.xlabel(r'$\log(\Sigma / \bar{\Sigma})$') + P.ylabel(r'$\mathcal{P}_\Sigma$') + if put_title: + P.title(title) + values, edges, _ = P.hist(dcoldens[mask_flat], + nb_bin_hist, range=(-1, 3), + weights = np.ones(nb_cells) / nb_cells) + centers = 0.5 * (edges[1:] + edges[:-1]) + + # Variance + var = np.var(dcoldens[mask_flat]) + # Compute the slope of the right part of the histogramm + mask_fit = (centers > 0.) & (centers < 1.25) & (values > 0) + if (np.sum(mask_fit > 0)): + (a, b, rho, _, stderr) = linregress(centers[mask_fit], np.log10(values[mask_fit])) + P.plot(centers, 10**(a*centers + b), '--', linewidth=2) + print("pdf a=%e, b=%e, rho=%e, var=%e"% (a, b, rho, var)) + try : + beta = int(tag.split('_')[1][4:]) + except ValueError : + beta = 0 + fit = {'beta': beta,'slope': a, 'origin': b, 'correlation': rho, 'stderr':stderr, 'var':var} + + return fit + +def plot_dcsdrho(fluct_maps, mask_flat, put_title=False, title='', nb_bin_hist=50, tag=''): + drho = fluct_maps['rho'].flatten() + dcs = fluct_maps['cs'].flatten() + + nb_points = np.sum(mask_flat) + size_hist = int(np.sqrt(nb_points)) + + P.hist2d(np.log10(drho[mask_flat]), np.log10(dcs[mask_flat]), (size_hist, size_hist), norm=mpl.colors.LogNorm()) + P.xlabel(r'$\log(\rho / \bar{\rho})$') + P.ylabel(r'$\log(c_s / \bar{c_s})$') + P.ylim([-0.6, 0.6]) + P.xlim([-1, 1.]) + + (a, b, rho, _, stderr) = linregress(np.log10(drho[mask_flat]), np.log10(dcs[mask_flat])) + P.plot(np.log10(drho[mask_flat]), a*np.log10(drho[mask_flat]) + b, '--', linewidth=2) + print("cs/rho a=%e, b=%e, rho=%e"% (a, b, rho)) + + try : + beta = int(tag.split('_')[1][4:]) + except ValueError : + beta = 0 + fit = {'beta': beta,'slope': a, 'origin': b, 'correlation': rho, 'stderr':stderr} + +def disk_pdf(path, num, maps_disk, + pos_star=[1., 1.], force=False, interactive=False, + nb_bin_hist=50, tag='', + rad_min=0.075, rad_max=0.3, put_title=True, + do_speed=True): # Load property file name_prop = path + "/prop_disk_" + str(num).zfill(5) + ".save" @@ -853,9 +906,11 @@ def disk_pdf( print("maps file loaded") # Properties - time = prop_disk["time"] - im_extent = maps_disk["im_extent"] - title = tag.split("_")[1] + " t=" + str(time)[0:5] + " (code)" + time = prop_disk['time'] + im_extent = maps_disk['im_extent'] + title = None + if (put_title): + title = tag.split('_')[1] + ' t='+ str(time)[0:5] +' (code)' # Load coldens coldens_map = maps_disk["coldens_z"] @@ -883,16 +938,14 @@ def disk_pdf( # Mask selecting the zone of interest mask_map = (rr > rad_min) & (rr < rad_max) - mask_flat = mask_map.flatten() # Additionnal maps rho_map = maps_disk["rho_z"] cs_map = np.sqrt(maps_disk["T_z"]) - vx_map = maps_disk["vx_z"] - vy_map = maps_disk["vy_z"] - v_map = np.sqrt(vx_map ** 2 + vy_map ** 2) - v_kepl = np.sqrt(1.0 / rr) + vx_map = maps_disk['vx_z'] + vy_map = maps_disk['vy_z'] + v_map = np.sqrt(vx_map**2 + vy_map**2) xx_star = xx - pos_star[0] yy_star = yy - pos_star[1] vrad_map = (vx_map * xx_star + vy_map * yy_star) / rr @@ -919,7 +972,7 @@ def disk_pdf( # save for later # prop_disk[cur_map + '_rmean'] = mean_bin - + # Add value for borders mean_bin = np.concatenate(([mean_bin[0]], mean_bin)) @@ -931,45 +984,29 @@ def disk_pdf( avg_maps[cur_map] = np.reshape(avg_flat, rr.shape) # Select zone of interest - avg_maps[cur_map][np.logical_not(mask_map)] = np.nan + # avg_maps[cur_map][np.logical_not(mask_map)] = np.nan # Compute fluctuation fluct_maps[cur_map] = map_arr / avg_maps[cur_map] + if not interactive: + plot_fluctuations_map(path, num, maps, fluct_maps, avg_maps, im_extent, mask_map, put_title, title, + nb_bin_hist, tag, prop_disk) + + # Save on disk + f = open(name_prop,'w') + pickle.dump(prop_disk, f) + f.close() + + return (xx, yy, fluct_maps) + +def plot_fluctuations_map(path, num, maps, fluct_maps, avg_maps, im_extent, mask_map, put_title, title, + nb_bin_hist, tag, prop_disk, interactive=False): + + mask_flat = mask_map.flatten() + # Sigma-PDF - dcoldens = np.log10(fluct_maps['coldens']).flatten() - - nb_cells = np.sum(mask_flat) - P.grid() - P.yscale('log') - P.ylim([0.5 / nb_cells, 1.]) - P.xlabel(r'$\log(\Sigma / \bar{\Sigma})$') - P.ylabel(r'$\mathcal{P}_\Sigma$') - if put_title: - P.title(title) - values, edges, _ = P.hist( - dcoldens[mask_flat], - nb_bin_hist, - range=(-1, 3), - weights=np.ones(nb_cells) / nb_cells, - ) - centers = 0.5 * (edges[1:] + edges[:-1]) - - # Variance - var = np.var(dcoldens[mask_flat]) - # Compute the slope of the right part of the histogramm - mask_fit = (centers > 0.) & (centers < 1.25) & (values > 0) - if (np.sum(mask_fit > 0)): - (a, b, rho, _, stderr) = linregress(centers[mask_fit], np.log10(values[mask_fit])) - P.plot(centers, 10**(a*centers + b), '--', linewidth=2) - print("a=%e, b=%e, rho=%e, var=%e"% (a, b, rho, var)) - try : - beta = int(tag.split('_')[1][4:]) - except ValueError : - beta = 0 - fit = {'beta': beta,'slope': a, 'origin': b, 'correlation': rho, 'stderr':stderr, 'var':var} - prop_disk['fit'] = fit - + fit = sigma_pdf(fluct_maps, mask_flat, put_title, title, nb_bin_hist, tag) if interactive: P.figure() else: @@ -999,49 +1036,25 @@ def disk_pdf( P.savefig(path + '/drho_hist_' + tag + '_' + str(num).zfill(5) + out_ext) P.close() - - # Derived quantities - dcs = fluct_maps['cs'].flatten() - dv = fluct_maps['v'].flatten() - dvaz_kepl = abs(maps['vaz'] - v_kepl) / v_kepl - fluct_maps['vaz_kepl'] = dvaz_kepl + dcoldens = np.log10(fluct_maps['coldens']).flatten() + dcs = fluct_maps['cs'].flatten() + dv = fluct_maps['v'].flatten() dmach = abs(maps['v'] - avg_maps['v']) / maps['cs'] dmach[dmach > 10.] = 10. fluct_maps['mach'] = dmach dmach_mean = np.mean(dmach[mask_map & (dmach < 5.)].flatten()) dcoldens_max = np.max(dcoldens[mask_flat]) - + # Save general properties prop_disk['dcoldens_max'] = dcoldens_max prop_disk['dmach_mean'] = dmach_mean print("dmach_mean = {}".format(dmach_mean)) - # Fluctuations plots # dcs = f(drho) - P.hist2d( - np.log10(drho[mask_flat]), - np.log10(dcs[mask_flat]), - (1000, 1000), - norm=mpl.colors.LogNorm(), - ) - P.xlabel(r"$\log(\rho / \bar{\rho})$") - P.ylabel(r"$\log(c_s / \bar{c_s})$") - P.ylim([-0.6, 0.6]) - P.xlim([-1, 1.]) - - (a, b, rho, _, stderr) = linregress(np.log10(drho[mask_flat]), np.log10(dcs[mask_flat])) - P.plot(np.log10(drho[mask_flat]), a*np.log10(drho[mask_flat]) + b, '--', linewidth=2) - print("cs/rho a=%e, b=%e, rho=%e"% (a, b, rho)) - - try : - beta = int(tag.split('_')[1][4:]) - except ValueError : - beta = 0 - fit = {'beta': beta,'slope': a, 'origin': b, 'correlation': rho, 'stderr':stderr} - prop_disk['fit_cs'] = fit + prop_disk['fit_cs'] = plot_dcsdrho(fluct_maps, mask_flat, put_title, title, nb_bin_hist, tag) P.colorbar() if interactive: @@ -1096,15 +1109,6 @@ def disk_pdf( vmin=-2., vmax=2.) label = r'$log(\rho/\bar{\rho})$' - elif cur_map == 'vaz_kepl': - im = P.imshow(fluct_map, - extent=im_extent, - origin='lower', - cmap='RdBu_r', - norm=mpl.colors.LogNorm(), - vmax=1., - vmin=0.01) - label = r'$|v_\varphi - v_{kepl}|/v_{kepl}$' elif cur_map == 'vaz': im = P.imshow(fluct_map, extent=im_extent, @@ -1135,15 +1139,9 @@ def disk_pdf( P.savefig(name_im) P.close() - # Save on disk - f = open(name_prop,'w') - pickle.dump(prop_disk, f) - f.close() - - def compare(path, runs, nums, path_out=None, force=False, interactive=False, Q_in_name=True, pdf=False, gamma_ad=5./3., tag=''): """ @@ -1168,7 +1166,7 @@ def compare(path, runs, nums, path_out=None, force=False, interactive=False, var_names = ['alpha_rey_mean', 'alpha_grav_mean', 'Q', 'time'] std_names = ['alpha_rey_mean', 'alpha_grav_mean', 'Q'] plot_names = ['alphaQ', 'alphaQall'] - + if Q_in_name: var_names.append('Q0') plot_names.append('alphaQ0') @@ -1187,20 +1185,20 @@ def compare(path, runs, nums, path_out=None, force=False, interactive=False, for name in std_names: std_data[name] = np.zeros(len(runs)) - + for i, run in enumerate(runs): path_run = path + "/" + run nb_outputs = 0 - + if(type(nums)==dict): nums_run = nums[run] else: nums_run = nums - + run_data = dict() for name in var_names: run_data[name] = [] - + for num in nums_run: try: # Load property file @@ -1227,11 +1225,11 @@ def compare(path, runs, nums, path_out=None, force=False, interactive=False, data = abs(prop_disk['alpha_grav_mean']) elif name == 'Q': data = prop_disk['Q_min'] - elif name == 'beta': + elif name == 'beta': data = fit['beta'] - elif name == 'kappa': + elif name == 'kappa': data = - fit['slope'] - elif name == 'var': + elif name == 'var': data = fit['var'] elif name == 'dmach': data = prop_disk['dmach_mean'] @@ -1242,10 +1240,10 @@ def compare(path, runs, nums, path_out=None, force=False, interactive=False, print(data) else: data = prop_disk[name] - + run_data[name].append(data) all_data[name].append(data) - + nb_outputs = nb_outputs + 1 print(run, num, nb_outputs) @@ -1272,7 +1270,7 @@ def compare(path, runs, nums, path_out=None, force=False, interactive=False, P.grid() print("Ploting "+ pname) - + if pname == 'alphaQ' or pname == 'alphaQ0': # alpha = f(Q) @@ -1284,11 +1282,11 @@ def compare(path, runs, nums, path_out=None, force=False, interactive=False, P.xlabel(r'$Q_{0}$') Q = var_data['Q0'] xerr = None - + P.yscale('log') P.ylim([1e-6,1.]) P.xlim([np.nanmin(Q)- 0.2, np.nanmax(Q) + 0.2]) - + P.errorbar(Q, abs(var_data['alpha_rey_mean']), xerr=xerr, yerr=std_data['alpha_rey_mean'], fmt='*--', label=r"$\bar{\alpha}_{Reynolds}$") # P.errorbar(Q, abs(var_data['alpha_grav_mean']), xerr=xerr, yerr=std_data['alpha_grav_mean'], @@ -1298,7 +1296,7 @@ def compare(path, runs, nums, path_out=None, force=False, interactive=False, P.ylabel(r'$\bar{\alpha}$') elif pname == 'kappabeta': P.errorbar(var_data['beta'], var_data['kappa'], yerr=std_data['kappa'], fmt='*') - + P.ylabel(r'$\kappa = - dlog\mathcal{P}_\Sigma / d\log\left(\Sigma/\overline{\Sigma}\right)$') P.xlabel(r'$\beta$') P.xlim([None, np.max(var_data['beta']) + 0.3]) @@ -1346,8 +1344,8 @@ def compare(path, runs, nums, path_out=None, force=False, interactive=False, P.xlabel(r'$<(v - \bar{v}) / c_s>$') P.ylabel(r'$\exp(Var(\log(\Sigma / \bar{\Sigma}))$') - if pname == 'alphaQall' : - Q = all_data['Q'] + if pname == 'alphaQall' : + Q = all_data['Q'] alpha_rey = abs(all_data['alpha_rey_mean']) time = all_data['time'] P.yscale('log') @@ -1451,12 +1449,12 @@ def evolution(path, nums, force=False, interactive=False, pdf=False, tag=''): except KeyError as e: print('[{}, {}] Key error {}'.format(path, num, e)) var_data[name][i] = np.nan - + time = var_data['time'] for pname in plot_names: - - if pname == 'alpha': + + if pname == 'alpha': P.yscale('log') P.plot(time, abs(var_data['alpha_rey_mean']), 'o-.', label=r"$\bar{\alpha}_{Reynolds}$") P.plot(time, abs(var_data['alpha_grav_mean']), '*--', label=r"$\bar{\alpha}_{grav}$") @@ -1465,7 +1463,7 @@ def evolution(path, nums, force=False, interactive=False, pdf=False, tag=''): P.plot(time, var_data['Q_min'], 'o-.', label=r"$Q_{min}$") P.plot(time, var_data['Q_mean'], '*--', label=r"$\overline{Q}$") P.legend() - elif pname == 'mass': + elif pname == 'mass': P.plot(time, var_data['mass_disk'], 'o-.', label=r"$M_{disk}$") P.plot(time, var_data['mass_box'], '*--', label=r"$M_{box}$") P.legend() @@ -1474,7 +1472,7 @@ def evolution(path, nums, force=False, interactive=False, pdf=False, tag=''): if pname =='dcoldens_max': P.plot(time, 1.5 * np.ones(len(time)), 'r-', lw=4) - + P.ylabel(plot_labels[pname]) P.xlabel(r'time (code)') P.grid() @@ -1484,7 +1482,7 @@ def evolution(path, nums, force=False, interactive=False, pdf=False, tag=''): P.tight_layout(pad=1) P.savefig(path + '/' + pname + 'time_' + tag + out_ext) P.close() - + def clump_study(path, num, path_out, tag, force=False, rho_thres=1e3, lvl_thres=10): name = tag + '_' + str(num).zfill(5) hop_save = name + '_hop' + '_prop_struct.save' @@ -1514,7 +1512,7 @@ def clump_study(path, num, path_out, tag, force=False, rho_thres=1e3, lvl_thres= print("Mpas files loaded") P.scatter(pos[:, 0],pos[:, 1], c = color) - + im = P.imshow(m['coldens_z'], extent=m['im_extent'], origin='lower', @@ -1541,3 +1539,142 @@ def clump_study(path, num, path_out, tag, force=False, rho_thres=1e3, lvl_thres= # P.imshow(np.log10(m['coldens_z']), extent=m['im_extent'], cmap='plasma') # P.savefig(path_out + '/clumps_size_' + tag + '_' + str(num).zfill(5) + out_ext) + + + +class interactive_plots: + + def onbuttonrelease(self, event): + """Deal with click events""" + button = ['left','middle','right'] + toolbar = P.get_current_fig_manager().toolbar + if toolbar.mode=='zoom rect' and event.inaxes == self.ax_col: + print("zooming ") + xlim = self.ax_col.get_xlim() + ylim = self.ax_col.get_ylim() + self.reset_mask() + elif self.add_mask and event.inaxes == self.ax_col: + self.plot_side() + P.draw() + + def onbuttonpress(self, event): + """Deal with click events""" + button = ['left','middle','right'] + toolbar = P.get_current_fig_manager().toolbar + if toolbar.mode!='': + print("You clicked on something, but toolbar is in mode {:s}.".format(toolbar.mode)) + print(self.add_mask) + if self.add_mask and toolbar.mode=='' and event.inaxes == self.ax_col: + ix, iy = event.xdata, event.ydata + print("Add patch {}, {}".format(ix, iy)) + xlim = self.ax_col.get_xlim() + ylim = self.ax_col.get_ylim() + radius = 0.05 * min(abs(xlim[1] - xlim[0]), abs(ylim[1] - ylim[0])) + circle = mpatches.Circle([ix, iy], radius, color='black', alpha=0.1, ec="none") + self.circles.append(circle) + self.ax_col.add_artist(circle) + self.ax_col.draw_artist(circle) + self.patch_mask = self.patch_mask | ((self.xx - ix)**2 + (self.yy -iy)**2 < radius**2) + #self.plot_side() + + def onkeypress(self, event): + """whenever a key is pressed""" + if not event.inaxes: + return + if event.key == 't': + self.add_mask = not self.add_mask + print("Add mode is {}".format(self.add_mask)) + elif event.key == 'r': + self.reset_mask() + + def plot_side(self): + if (self.add_mask): + mask = (self.patch_mask & self.mask).flatten() + else: + mask = self.mask.flatten() + self.ax_gamma.clear() + P.sca(self.ax_gamma) + plot_dcsdrho(self.fluct_maps, mask, tag=self.tag) + + self.ax_pdf.clear() + P.sca(self.ax_pdf) + sigma_pdf(self.fluct_maps, mask, tag=self.tag, nb_bin_hist=self.args.pdf_nb_bin) + + def reset_mask(self): + xlim = self.ax_col.get_xlim() + ylim = self.ax_col.get_ylim() + self.mask = (self.xx >= xlim[0]) & (self.xx <= xlim[1]) & (self.yy >= ylim[0]) & (self.yy <= ylim[1]) + self.patch_mask = np.full(self.mask.shape, False) + for circle in self.circles: + circle.remove() + self.circles = [] + self.plot_side() + P.draw() + + + + def __init__(self, args, path, num, + maps_disk=None, + tag='', + set_lim=True) : + """ + Interactive plotting + + Parameters + ---------- + num output number + path path of the pipeline output + + """ + self.args = args + self.add_mask = False + self.circles = [] + self.tag = tag + + path_out = path + + # Load maps file + print("load maps file") + name_maps = path + '/maps_disk' + '_' + tag + '_' + format(num,'05') + '.save' + + if maps_disk is None: + if (len(glob.glob(name_maps)) == 0): + raise IOError('no pickle file for disk maps {}. Run make_image_disk() first'.format(name_maps)) + f = open(name_maps,'r') + maps_disk = pickle.load(f) + f.close() + print("maps file loaded") + + im_extent = maps_disk['im_extent'] + + fig = P.figure(); + self.ax_col = P.subplot(1, 2, 1) + coldens = maps_disk['coldens_z'] + im = self.ax_col.imshow(coldens, + extent=im_extent, + origin='lower', + norm=mpl.colors.LogNorm()) + if set_lim: + im.set_clim(0.01, 100) + self.ax_col.set_xlabel(r'$x$') + self.ax_col.set_ylabel(r'$y$') + + self.xx, self.yy, self.fluct_maps = disk_pdf(path, num, maps_disk, tag=self.tag, force=True, put_title=False, interactive=True) + coord_flat = zip(self.xx.flatten(), self.yy.flatten()) + + self.ax_gamma = P.subplot(2, 2, 2) + self.ax_pdf = P.subplot(2,2,4) + + xlim = self.ax_col.get_xlim() + ylim = self.ax_col.get_ylim() + self.mask = (self.xx >= xlim[0]) & (self.xx <= xlim[1]) & (self.yy >= ylim[0]) & (self.yy <= ylim[1]) + self.patch_mask = np.full(self.mask.shape, False) + + self.plot_side() + + fig.canvas.mpl_connect('button_release_event', self.onbuttonrelease) + fig.canvas.mpl_connect('button_press_event', self.onbuttonpress) + fig.canvas.mpl_connect('key_press_event', self.onkeypress) + + P.tight_layout() + P.show() diff --git a/pipeline_disk.py b/pipeline_disk.py index a2e682f..c41cd12 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -109,26 +109,23 @@ parser.add_argument("--print_outputs", help="print names of outputs", action='store_true') -parser.add_argument( - "--fft", help="use quick and dirty fft rendering", action="store_true" -) -parser.add_argument( - "--images", - nargs="*", - default=["coldens", "rho", "speed", "Q", "T"], - choices=["coldens", "rho", "speed", "Q", "T", "levels", "cpu"], -) -parser.add_argument( - "--axes", nargs="*", default=["x", "y", "z"], choices=["x", "y", "z"] -) -parser.add_argument("--zoom", help="zoom", type=float, default=2.0) -parser.add_argument( - "-ms", - "--mapsize", - help="size of the maps created in he map mode (in pixel)", - type=int, - default=1024, -) +parser.add_argument("--fft", + help="use quick and dirty fft rendering", + action='store_true') +parser.add_argument("--images", nargs='*', + default=['coldens', 'rho', 'speed', 'Q', 'T'], + choices=['coldens', 'rho', 'speed', 'Q', 'T', 'levels', 'cpu', 'jeans', 'jeans_ratio']) +parser.add_argument("--axes", nargs='*', + default=['x', 'y', 'z'], + choices=['x', 'y', 'z']) +parser.add_argument("--zoom", + help="zoom", + type=float, + default=2.) +parser.add_argument("-ms", "--mapsize", + help="size of the maps created in he map mode (in pixel)", + type=int, + default=1024) parser.add_argument("--nb_bin", @@ -206,8 +203,8 @@ if args.format == 'pdf': if args.beamer: - dp.P.rcParams['font.family'] = 'sans-serif' - dp.P.rcParams['figure.figsize'] = (7, 4.5) + dp.P.rcParams['font.family'] = 'sans-serif' + dp.P.rcParams['figure.figsize'] = (7, 4.5) # Plot properties dp.P.rcParams['image.cmap'] = args.colormap @@ -221,148 +218,166 @@ me.P = dp.P # List of id that were successfully computed run_succeded = {} + +# Care for dependencies +images = args.images +if 'jeans_ratio' in images and not 'jeans' in images : + images = ['jeans'] + images +if 'jeans_ratio' in images and not 'levels' in images : + images.append('levels') +if 'jeans' in images and not 'T' in images : + images.append('T') +if ('T' in images or 'jeans' in images) and not 'rho' in images : + images.append('rho') +if 'Q' in images and not 'coldens' in images : + images.append('coldens') + +# Go through all runs for run in runs: - path_suffix = project + "/" + run - path_in = storage_in + path_suffix - path_out = storage_out + path_suffix + path_suffix = project + '/' + run + path_in = storage_in + path_suffix + path_out = storage_out + path_suffix - if args.tag == '': - tag_run = run - else: - tag_run = run + '_' + args.tag + if args.tag == '': + tag_run = run + else: + tag_run = run + '_' + args.tag - if not os.path.exists(path_out): - os.makedirs(path_out) - try: - copy(path_in + "/disk.nml", path_out) - copy(path_in + "/output_00001/compilation.txt", path_out) - except: - pass + if not os.path.exists(path_out): + os.makedirs(path_out) + try: + copy(path_in + '/disk.nml', path_out) + copy(path_in + '/output_00001/compilation.txt', path_out) + except: + pass - run_succeded[run] = [] + run_succeded[run] = [] - if args.which_outputs in ['all', 'time'] : - names = glob.glob(path_in + "/output_[0-9][0-9][0-9][0-9][0-9]") - nums_all = [int(n.split('/')[-1].split('_')[1]) for n in names] - nums_all = np.sort(nums_all) - if args.which_outputs == 'all': - nums = nums_all - else: - time = [dp.get_time(path_in, n) for n in nums_all] - nums = [n for i,n in enumerate(nums_all) if time[i] >= args.time_begin and time[i] < args.time_end] - else: - nums = range(first, last + 1, step) + if args.which_outputs in ['all', 'time'] : + names = glob.glob(path_in + "/output_[0-9][0-9][0-9][0-9][0-9]") + nums_all = [int(n.split('/')[-1].split('_')[1]) for n in names] + nums_all = np.sort(nums_all) + if args.which_outputs == 'all': + nums = nums_all + else: + time = [dp.get_time(path_in, n) for n in nums_all] + nums = [n for i,n in enumerate(nums_all) if time[i] >= args.time_begin and time[i] < args.time_end] + else: + nums = range(first, last + 1, step) - for num in nums: - failures = 0 - success = False + for num in nums: + failures = 0 + success = False - while not success: - try: - maps_disk = None - if args.print_outputs: - print("[{}, {}]".format(run, num)) - if args.maps: - print("[{}, {}] computing maps".format(run, num)) - maps_disk = dp.compute_image_data(path_in, num, - radius=rad, - path_out=path_out, - tag=tag_run, - map_size=args.mapsize, - force=args.force_redo, - axes_los=args.axes, - images=args.images, - pos_star=np.array([args.x, args.y, args.z]), - fft=args.fft) - print("[{}, {}] maps computed".format(run, num)) - if args.plot_maps: - print("[{}, {}] plotting maps".format(run, num)) - maps_disk = dp.plot_maps(path_out, num, - maps_disk=maps_disk, - axes_los=args.axes, - images=args.images, - tag=tag_run, - force=args.force_redo, - interactive=args.interactive, - put_title=(not args.beamer)) - print("[{}, {}] maps plotted".format(run, num)) - if args.disk: - print("[{}, {}] computing disk props".format(run, num)) - dp.disk_prop(path_in, num, path_out=path_out, - nb_bin=args.nb_bin, - binning=args.binning, - rad_ext=args.rad_ext, - force=args.force_redo, - pos_star=np.array([args.x, args.y, args.z])) - print("[{}, {}] disk_props computed".format(run, num)) - if args.plot_disk: - print("[{}, {}] plotting disk props".format(run, num)) - dp.plot_disk_prop(path_out, num, tag=tag_run, + while not success: + try: + maps_disk = None + if args.print_outputs: + print("[{}, {}]".format(run, num)) + if args.maps: + print("[{}, {}] computing maps".format(run, num)) + maps_disk = dp.compute_image_data(path_in, num, + radius=rad, + path_out=path_out, + tag=tag_run, + map_size=args.mapsize, + force=args.force_redo, + axes_los=args.axes, + images=images, + pos_star=np.array([args.x, args.y, args.z]), + fft=args.fft) + print("[{}, {}] maps computed".format(run, num)) + if args.plot_maps: + print("[{}, {}] plotting maps".format(run, num)) + maps_disk = dp.plot_maps(path_out, num, + maps_disk=maps_disk, + axes_los=args.axes, + images=images, + tag=tag_run, + force=args.force_redo, + interactive=args.interactive, + put_title=(not args.beamer)) + print("[{}, {}] maps plotted".format(run, num)) + if args.disk: + print("[{}, {}] computing disk props".format(run, num)) + dp.disk_prop(path_in, num, path_out=path_out, + nb_bin=args.nb_bin, + binning=args.binning, + rad_ext=args.rad_ext, force=args.force_redo, - interactive=args.interactive, - put_title=(not args.beamer)) - print("[{}, {}] disk props plotted".format(run, num)) - if args.pdf: - print("[{}, {}] computing pdf #1".format(run, num)) - dp.disk_pdf(path_out, num, maps_disk, - pos_star=np.array([args.x, args.y, args.z]), - force=args.force_redo, - tag=tag_run, - nb_bin_hist=args.pdf_nb_bin, - interactive=args.interactive, - put_title=(not args.beamer)) - print("[{}, {}] pdf #1 computed".format(run, num)) - if args.pdf2: - print("[{}, {}] computing pdf #2".format(run, num)) - me.get_pdf(path_in, num, path_out=path_out, - force=args.force_redo, - tag=tag_run) - print("[{}, {}] pdf #2 computed".format(run, num)) - if args.clump: - print("[{}, {}] computing clumps".format(run, num)) - dp.clump_study(path_in, num, path_out, tag_run) - print("[{}, {}] clumps computed".format(run, num)) - # If we are here, success ! - success = True - run_succeded[run].append(num) - except (ValueError, IOError, KeyError) as e: - print(e) - if args.watch and failures < args.allowed_failures: - failures = failures + 1 - print("Unable to proceed for run {} output {}. Trying again in {} s ({} tries remaining)".format(run, num, args.waiting_time, args.allowed_failures - failures)) - time.sleep(args.waiting_time) - elif args.skip: - break - else: - raise - if args.evolution: - print("[{}] plotting evolution".format(run)) - dp.evolution(path_out, run_succeded[run], - force=args.force_redo, - tag=args.tag, - interactive=args.interactive, - pdf=args.pdf or args.cpdf) - print("[{}] evolution plotted".format(run)) + pos_star=np.array([args.x, args.y, args.z])) + print("[{}, {}] disk_props computed".format(run, num)) + if args.plot_disk: + print("[{}, {}] plotting disk props".format(run, num)) + dp.plot_disk_prop(path_out, num, tag=tag_run, + force=args.force_redo, + interactive=args.interactive, + put_title=(not args.beamer)) + print("[{}, {}] disk props plotted".format(run, num)) + if args.pdf: + print("[{}, {}] computing pdf #1".format(run, num)) + dp.disk_pdf(path_out, num, maps_disk, + pos_star=np.array([args.x, args.y, args.z]), + force=args.force_redo, + tag=tag_run, + nb_bin_hist=args.pdf_nb_bin, + interactive=args.interactive, + put_title=(not args.beamer)) + print("[{}, {}] pdf #1 computed".format(run, num)) + if args.pdf2: + print("[{}, {}] computing pdf #2".format(run, num)) + me.get_pdf(path_in, num, path_out=path_out, + force=args.force_redo, + tag=tag_run) + print("[{}, {}] pdf #2 computed".format(run, num)) + if args.clump: + print("[{}, {}] computing clumps".format(run, num)) + dp.clump_study(path_in, num, path_out, tag_run) + print("[{}, {}] clumps computed".format(run, num)) + if args.interactive: + print("[{}, {}] Interactive session".format(run, num)) + maps_disk = dp.interactive_plots(args, path_out, num, tag=tag_run) + # If we are here, success ! + success = True + run_succeded[run].append(num) + except (ValueError, IOError, KeyError) as e: + print(e) + if(args.watch and failures < args.allowed_failures): + failures = failures + 1 + print("Unable to proceed for run {} output {}. Trying again in {} s ({} tries remaining)".format(run, num, args.waiting_time, args.allowed_failures - failures)) + time.sleep(args.waiting_time) + elif args.skip: + break + else: + raise + if args.evolution: + print("[{}] plotting evolution".format(run)) + dp.evolution(path_out, run_succeded[run], + force=args.force_redo, + tag=args.tag, + interactive=args.interactive, + pdf=args.pdf or args.cpdf) + print("[{}] evolution plotted".format(run)) if args.compare: - path_suffix = project - path = storage_out + path_suffix + path_suffix = project + path = storage_out + path_suffix - if (args.time_avg): - dp.compare(path, runs, run_succeded, - path_out=path + '/comp', - force=args.force_redo, - interactive=args.interactive, - Q_in_name=args.Q_in_name, - tag=args.tag, - pdf=args.pdf or args.cpdf) - else: - for i in range(first, last + 1, step): - # Select usable runs - runs_ok = [run for run in runs if i in run_succeded[run]] + if (args.time_avg): + dp.compare(path, runs, run_succeded, + path_out=path + '/comp', + force=args.force_redo, + interactive=args.interactive, + Q_in_name=args.Q_in_name, + tag=args.tag, + pdf=args.pdf or args.cpdf) + else: + for i in range(first, last + 1, step): + # Select usable runs + runs_ok = [run for run in runs if i in run_succeded[run]] - try: - if len(runs_ok) > 1: + try: + if len(runs_ok) > 1: dp.compare(path, runs_ok, i, path_out=path + '/comp', force=args.force_redo, @@ -370,9 +385,9 @@ if args.compare: Q_in_name=args.Q_in_name, tag=args.tag, pdf=args.pdf or args.cpdf) - except (KeyError, IOError) as e: - print(e) - if args.skip: + except (KeyError, IOError) as e: + print(e) + if (args.skip): pass - else: + else: raise