linear interpolation + regression

disk_postprocess.py

@@ -19,7 +19,7 @@ except:
     print("cPickle not found, using pickle")
     import pickle
 import tables
-
+from scipy.stats import linregress
 from pymses.sources.ramses import output
 from pymses.analysis import Camera, raytracing, slicing, splatting
 from pymses.filters import CellsToPoints
@@ -872,8 +872,9 @@ def disk_pdf(
     interactive=False,
     nb_bin_hist=50,
     tag="",
-    rad_min=0.08,
+    rad_min=0.05,
     rad_max=0.4,
+    put_title=True,
 ):
 
     # Load property file
@@ -900,7 +901,7 @@ def disk_pdf(
     print("maps file loaded")
 
     time = prop_disk["time"]
-    title = "t=" + str(time)[0:5] + " (code)"
+    title = tag.split("_")[1] + " t=" + str(time)[0:5] + " (code)"
 
     coldens_map = maps_disk["coldens_z"]
     im_extent = maps_disk["im_extent"]
@@ -945,39 +946,58 @@ def disk_pdf(
 
     # Histogramm : density fluctuation distribution function
     dcoldens = np.log10(coldens_map.flatten() / coldens_mean)
+    mask_hist = (rr_flat > rad_min) & (rr_flat < rad_max)
+    nb_cells = np.sum(mask_hist)
     P.grid()
     P.yscale("log")
+    P.ylim([0.5 / nb_cells, 1.0])
     P.xlabel(r"$\log(N / \bar{N})$")
     P.ylabel("Probability density")
-    P.title(title)
+    if put_title:
-    P.hist(
+        P.title(title)
-        dcoldens[(rr_flat > rad_min) & (rr_flat < rad_max)], nb_bin_hist, range=(-1, 3)
+    values, edges, _ = P.hist(
+        dcoldens[mask_hist],
+        nb_bin_hist,
+        range=(-1, 3),
+        weights=np.ones(nb_cells) / nb_cells,
     )
+    centers = 0.5 * (edges[1:] + edges[:-1])
+
+    # 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 interactive:
         pass
     else:
-        P.savefig(path + "/dcol_hist_" + str(num).zfill(5) + out_ext)
+        P.savefig(path + "/dcol_hist_" + tag + "_" + str(num).zfill(5) + out_ext)
         P.close()
 
     # Fluctuation map
     dcoldens_map = np.reshape(dcoldens, coldens_map.shape)
 
-    # cont = P.contour(coldens_mean_map,
-    #                  extent=im_extent,
-    #                  origin='lower',
-    #                  colors='k',
-    #                  linewidths=0.1)
-
     im = P.imshow(dcoldens_map, extent=im_extent, origin="lower", cmap="viridis")
 
-    P.title(title)
+    if put_title:
+        P.title(title)
     P.xlabel(r"$x$")
     P.ylabel(r"$y$")
 
     cbar = P.colorbar(im)
     cbar.set_label(r"$log(N/\bar{N})$ (code)")
-    name = path + "/dcoldensmap_" + format(num, "05")
+    name = path + "/dcoldensmap_" + tag + "_" + format(num, "05")
     name_im = name + out_ext
 
     if interactive:
@@ -987,7 +1007,7 @@ def disk_pdf(
         P.close()
 
 
-def compare(path, runs, num, force=False, interactive=False):
+def compare(path, runs, num, force=False, interactive=False, Q_in_name=True, pdf=False):
     """
     Compare properties of a disk in several simulations
 
@@ -1003,7 +1023,11 @@ def compare(path, runs, num, force=False, interactive=False):
     alpha_rey = np.zeros(len(runs))
     alpha_grav = np.zeros(len(runs))
     Q = np.zeros(len(runs))
-    Q0 = np.zeros(len(runs))
+    if Q_in_name:
+        Q0 = np.zeros(len(runs))
+    if pdf:
+        beta = np.zeros(len(runs))
+        slope = np.zeros(len(runs))
 
     for i, run in enumerate(runs):
         path_run = path + "/" + run
@@ -1021,12 +1045,17 @@ def compare(path, runs, num, force=False, interactive=False):
         alpha_rey[i] = prop_disk["alpha_rey_mean"]
         alpha_grav[i] = prop_disk["alpha_grav_mean"]
         Q[i] = prop_disk["Q_min"]
-        Q0[i] = float(run.split("_")[2][1:])
+        if Q_in_name:
+            Q0[i] = float(run.split("_")[2][1:])
+        if pdf:
+            fit = prop_disk["fit"]
+            beta[i] = fit["beta"]
+            slope[i] = fit["slope"]
 
     # Check if the output file exists, and exit if it is the case
     name_save = path + "/alphaQ_" + str(num).zfill(5) + out_ext
-    if not force and len(glob.glob(name_save)) != 0:
+    # if (not force and len(glob.glob(name_save)) !=0):
-        return
+    #     return
 
     title = "t=" + str(time[0]) + " (code)"
 
@@ -1048,23 +1077,57 @@ def compare(path, runs, num, force=False, interactive=False):
         P.savefig(path + "/alphaQ_" + str(num).zfill(5) + out_ext)
         P.close()
 
-    # alpha = f(Q0)
+    if Q_in_name:
-    P.yscale("log")
+        # alpha = f(Q0)
-    P.ylim([1e-7, 1.0])
+        P.yscale("log")
-    P.grid()
+        P.ylim([1e-7, 1.0])
+        P.grid()
 
-    P.plot(Q0, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$")
+        P.plot(Q0, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$")
-    P.plot(Q0, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$")
+        P.plot(Q0, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$")
 
-    P.legend()
+        P.legend()
-    P.ylabel(r"$\bar{\alpha}$")
+        P.ylabel(r"$\bar{\alpha}$")
-    P.xlabel(r"$Q_{0}$")
+        P.xlabel(r"$Q_{0}$")
 
-    if interactive:
+        if interactive:
-        P.figure()
+            P.figure()
-    else:
+        else:
-        P.savefig(path + "/alphaQ0_" + str(num).zfill(5) + out_ext)
+            P.savefig(path + "/alphaQ0_" + str(num).zfill(5) + out_ext)
-        P.close()
+            P.close()
+
+    if pdf:
+        # slope of the pdf = f(beta)
+        P.grid()
+        P.plot(beta, slope, "o-.")
+
+        P.legend()
+        P.ylabel(r"$d\log\% / d\logN$")
+        P.xlabel(r"$\beta$")
+
+        if interactive:
+            P.figure()
+        else:
+            P.savefig(path + "/dcolslopebeta_" + str(num).zfill(5) + out_ext)
+            P.close()
+
+        # alpha = f(beta)
+        P.yscale("log")
+        P.ylim([1e-7, 1.0])
+        P.grid()
+
+        P.plot(beta, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$")
+        P.plot(beta, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$")
+
+        P.legend()
+        P.ylabel(r"$\bar{\alpha}$")
+        P.xlabel(r"$\beta_{cool}$")
+
+        if interactive:
+            P.figure()
+        else:
+            P.savefig(path + "/alphabeta_" + str(num).zfill(5) + out_ext)
+            P.close()
 
 
 def evolution(path, nums, force=False, interactive=False):

pipeline_disk.py

@@ -224,4 +224,12 @@ if args.compare:
     path = storage_out + path_suffix
 
     for i in range(first, last + 1, step):
-        dp.compare(path, runs, i, force=args.force_redo, interactive=args.interactive)
+        dp.compare(
+            path,
+            runs,
+            i,
+            force=args.force_redo,
+            interactive=args.interactive,
+            Q_in_name=(not args.pdf),
+            pdf=args.pdf,
+        )