From aa8ba76162071f819cf4754b948f0ee4eadf256a Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Mon, 20 May 2019 10:56:31 +0200 Subject: [PATCH] Time average, other improvments --- disk_postprocess.py | 124 +++++++++++++++++++++++++++----------------- pipeline_disk.py | 34 +++++++----- 2 files changed, 99 insertions(+), 59 deletions(-) diff --git a/disk_postprocess.py b/disk_postprocess.py index e151076..27f2085 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -152,8 +152,8 @@ def make_image_aux( directory = path # Checking for existing files - name = directory + "/coldens_z" + "_" + tag + "_" + format(num, "05") + out_ext - if len(glob.glob(name)) == 1 and not force: + name_save = directory + "/maps_disk" + "_" + tag + "_" + format(num, "05") + ".save" + if len(glob.glob(name_save)) == 1 and not force: return # Prepare saving data @@ -453,9 +453,6 @@ def make_image_aux( P.close() if save_data: - name_save = ( - directory + "/maps_disk" + "_" + tag + "_" + format(num, "05") + ".save" - ) f = open(name_save, "w") pickle.dump(maps_disk, f) f.close() @@ -961,7 +958,7 @@ def disk_pdf( P.yscale("log") P.ylim([0.5 / nb_cells, 1.0]) P.xlabel(r"$\log(N / \bar{N})$") - P.ylabel("Probability density") + P.ylabel(r"\mathcal{P}_N") if put_title: P.title(title) values, edges, _ = P.hist( @@ -974,19 +971,20 @@ def disk_pdf( # Compute the slope of the right part of the histogramm mask_fit = (centers > 0.25) & (centers < 1.25) & (values > 0) - (a, b, rho, _, _) = linregress(centers[mask_fit], np.log10(values[mask_fit])) - print("a=%e, b=%e, rho=%e" % (a, b, rho)) - fit = { - "beta": int(tag.split("_")[1][4:]), - "slope": a, - "origin": b, - "correlation": rho, - } - P.plot(centers, 10 ** (a * centers + b), "--") - f = open(name_prop, "w") - prop_disk["fit"] = fit - pickle.dump(prop_disk, f) - f.close() + if np.sum(mask_fit > 0): + (a, b, rho, _, _) = linregress(centers[mask_fit], np.log10(values[mask_fit])) + P.plot(centers, 10 ** (a * centers + b), "--") + print("a=%e, b=%e, rho=%e" % (a, b, rho)) + fit = { + "beta": int(tag.split("_")[1][4:]), + "slope": a, + "origin": b, + "correlation": rho, + } + f = open(name_prop, "w") + prop_disk["fit"] = fit + pickle.dump(prop_disk, f) + f.close() if interactive: pass @@ -1020,6 +1018,7 @@ def compare( path, runs, nums, + path_out=None, force=False, interactive=False, Q_in_name=True, @@ -1039,7 +1038,7 @@ def compare( if type(nums) == int: nums = [nums] - nums_name = "_".join(str(num).zfill(5) for num in nums) + nums_name = "_".join(str(num).zfill(5) for num in [nums[0], nums[-1]]) # Initialize arrays alpha_rey = np.zeros(len(runs)) @@ -1053,25 +1052,49 @@ def compare( for i, run in enumerate(runs): path_run = path + "/" + run + nb_outputs = 0 - for j, num in enumerate(nums): - # Load property file - name_save = path_run + "/prop_disk_" + str(num).zfill(5) + ".save" + for num in nums: + try: + # Load property file + name_save = path_run + "/prop_disk_" + str(num).zfill(5) + ".save" - # Check if the properties file exists - if len(glob.glob(name_save)) == 0: - raise ("no pickle file for disk properties. Run disk_prop() first") - f = open(name_save, "r") - prop_disk = pickle.load(f) - f.close() + # Check if the properties file exists + if len(glob.glob(name_save)) == 0: + raise IOError( + "no pickle file for disk properties {}. Run disk_prop() first".format( + name_save + ) + ) + f = open(name_save, "r") + prop_disk = pickle.load(f) + f.close() - alpha_rey[i] = alpha_rey[i] + prop_disk["alpha_rey_mean"] / len(nums) - alpha_grav[i] = alpha_grav[i] + prop_disk["alpha_grav_mean"] / len(nums) - Q[i] = Q[i] + prop_disk["Q_min"] / len(nums) + alpha_rey[i] = alpha_rey[i] + abs(prop_disk["alpha_rey_mean"]) + alpha_grav[i] = alpha_grav[i] + abs(prop_disk["alpha_grav_mean"]) + Q[i] = Q[i] + prop_disk["Q_min"] + if pdf: + fit = prop_disk["fit"] + beta[i] = fit["beta"] + slope[i] = slope[i] + fit["slope"] + + nb_outputs = nb_outputs + 1 + print(run, num, nb_outputs) + + except (IOError, KeyError) as e: + print(run, num, e, nb_outputs) + + if nb_outputs > 0: + alpha_rey[i] = alpha_rey[i] / nb_outputs + alpha_grav[i] = alpha_grav[i] / nb_outputs + Q[i] = Q[i] / nb_outputs if pdf: - fit = prop_disk["fit"] - beta[i] = fit["beta"] - slope[i] = slope[i] + fit["slope"] / len(nums) + slope[i] = slope[i] / nb_outputs + else: + for array in [alpha_rey, alpha_grav, Q]: + array[i] = np.nan + if pdf: + slope[i] = np.nan if Q_in_name: Q0[i] = float(run.split("_")[2][1:]) @@ -1096,7 +1119,7 @@ def compare( if interactive: P.figure() else: - P.savefig(path + "/alphaQ_" + nums_name + out_ext) + P.savefig(path_out + "/alphaQ_" + nums_name + out_ext) P.close() if Q_in_name: @@ -1115,7 +1138,7 @@ def compare( if interactive: P.figure() else: - P.savefig(path + "/alphaQ0_" + nums_name + out_ext) + P.savefig(path_out + "/alphaQ0_" + nums_name + out_ext) P.close() if pdf: @@ -1130,7 +1153,7 @@ def compare( if interactive: P.figure() else: - P.savefig(path + "/dcolslopebeta_" + nums_name + out_ext) + P.savefig(path_out + "/dcolslopebeta_" + nums_name + out_ext) P.close() # alpha = f(beta) @@ -1152,7 +1175,7 @@ def compare( if interactive: P.figure() else: - P.savefig(path + "/alphabeta_" + nums_name + out_ext) + P.savefig(path_out + "/alphabeta_" + nums_name + out_ext) P.close() @@ -1190,14 +1213,20 @@ def evolution(path, nums, force=False, interactive=False, pdf=False): f.close() time[i] = prop_disk["time"] - alpha_rey[i] = prop_disk["alpha_rey_mean"] - alpha_grav[i] = prop_disk["alpha_grav_mean"] - Qmin[i] = prop_disk["Q_min"] - Qmean[i] = prop_disk["Q_mean"] - mass_disk[i] = prop_disk["mass_disk"] - mass_box[i] = prop_disk["mass_box"] + + try: + alpha_rey[i] = prop_disk["alpha_rey_mean"] + alpha_grav[i] = prop_disk["alpha_grav_mean"] + Qmin[i] = prop_disk["Q_min"] + Qmean[i] = prop_disk["Q_mean"] + mass_disk[i] = prop_disk["mass_disk"] + mass_box[i] = prop_disk["mass_box"] + except: + for array in [alpha_rey, alpha_grav, Qmin, Qmean, mass_disk, mass_box]: + array[i] = np.nan + if pdf: - slope = prop_disk["fit"]["slope"] + slope[i] = prop_disk["fit"]["slope"] # Check if the output file exists, and exit if it is the case name_save = path + "/alpha_time" + out_ext @@ -1255,9 +1284,10 @@ def evolution(path, nums, force=False, interactive=False, pdf=False): # PDF if pdf: P.grid() + print(time, slope) P.plot(time, slope, "o-.") P.legend() - P.ylabel(r"$d\log\% / d\logN$") + P.ylabel(r"$d\log\mathcal{P}_{N} / d\logN$") P.xlabel(r"time (code)") if interactive: P.figure() diff --git a/pipeline_disk.py b/pipeline_disk.py index 3ee93f2..704c638 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -238,14 +238,15 @@ if args.compare: path = storage_out + path_suffix if args.time_avg: - # Select output availables for all runs - output_ok = reduce(np.intersect1d, [run_succeded[run] for run in runs]) + # Select output availables for at least on run + output_ok = reduce(np.union1d, [run_succeded[run] for run in runs]).astype(int) print(output_ok) if len(output_ok) >= 1: dp.compare( path, runs, output_ok, + path_out=path + "/comp", force=args.force_redo, interactive=args.interactive, Q_in_name=(not args.pdf), @@ -255,13 +256,22 @@ if args.compare: 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 len(runs_ok) > 1: - dp.compare( - path, - runs_ok, - i, - force=args.force_redo, - interactive=args.interactive, - Q_in_name=(not args.pdf), - pdf=args.pdf, - ) + + try: + if len(runs_ok) > 1: + dp.compare( + path, + runs_ok, + i, + path_out=path + "/comp", + force=args.force_redo, + interactive=args.interactive, + Q_in_name=(not args.pdf), + pdf=args.pdf, + ) + except (KeyError, IOError) as e: + print(e) + if args.skip: + pass + else: + raise