rework time-series plotting; plot results for multiple diffusion constants

Author: jonaspleyer

cr_bayesian_optim/fractal_dim.py

@@ -123,45 +123,73 @@ def produce_options():
 def fractal_dim_over_time():
     options = produce_options()
 
-    cells, _ = load_or_compute_full(options)
-
-    iterations = sorted(cells.keys())[::4]
-    colony_diam = []
-    dims_mean = []
-    dims_std = []
-    for i in tqdm(iterations, desc="Calculating dim(t)"):
-        pos = np.array([c[0].mechanics.pos for c in cells[i].values()])
-
-        # Calculate Diameter of colony with convex hull
-        hull = sp.spatial.ConvexHull(pos)
-        hull_points = pos[hull.vertices]
-
-        diam = 0
-        n_points = int(np.ceil(len(hull_points) / 2))
-        for p in hull_points[:n_points]:
-            d = np.linalg.norm(hull_points - p, axis=1)
-            diam = max(diam, np.max(d))
-        colony_diam.append(diam)
-
-        _, _, popt, pcov = calculate_fractal_dim_for_pos(pos, options, None)
-        dims_mean.append(-popt[0])
-        dims_std.append(pcov[0, 0] ** 0.5)
-
-    t = np.array(iterations) * options.time.dt / 60
-    y1 = np.array(dims_mean)
-    y1_err = np.array(dims_std)
-    y2 = np.array(colony_diam) / 1000
+    t = []
+    y1 = []
+    y1_err = []
+    y2 = []
+
+    diffusion_constants = [80, 5, 0.5]
+    for diffusion_constant in diffusion_constants:
+        options.domain.diffusion_constant = diffusion_constant
+        cells, _ = load_or_compute_full(options)
+
+        iterations = sorted(cells.keys())[::4]
+        colony_diam = []
+        dims_mean = []
+        dims_std = []
+        for i in tqdm(iterations, desc="Calculating dim(t)"):
+            pos = np.array([c[0].mechanics.pos for c in cells[i].values()])
+
+            # Calculate Diameter of colony with convex hull
+            hull = sp.spatial.ConvexHull(pos)
+            hull_points = pos[hull.vertices]
+
+            diam = 0
+            n_points = int(np.ceil(len(hull_points) / 2))
+            for p in hull_points[:n_points]:
+                d = np.linalg.norm(hull_points - p, axis=1)
+                diam = max(diam, np.max(d))
+            colony_diam.append(diam)
+
+            _, _, popt, pcov = calculate_fractal_dim_for_pos(pos, options, None)
+            dims_mean.append(-popt[0])
+            dims_std.append(pcov[0, 0] ** 0.5)
+
+        t.append(np.array(iterations) * options.time.dt / 60)
+        y1.append(np.array(dims_mean))
+        y1_err.append(np.array(dims_std))
+        y2.append(np.array(colony_diam) / 1000)
 
     fig, ax = plt.subplots(figsize=(8, 8))
 
     # Introduce new y-axis for colony size
     ax2 = ax.twinx()
     ax2.set_ylabel("Diameter [mm]")
-    ax2.plot(t, y2, label="Colony Size", linestyle="--", color=COLOR5)
+    ax2.set_yscale("log")
+    for i in range(len(t)):
+        ax2.plot(
+            t[i], y2[i], label="Colony Size", linestyle="--", color=COLOR5, linewidth=2
+        )
 
     # Plot Fractal Dimension
-    ax.plot(t, y1, label="dim", color=COLOR1)
-    ax.fill_between(t, y1 - y1_err, y1 + y1_err, color=COLOR3, alpha=0.3)
+    for i in range(len(t)):
+        ax.plot(t[i], y1[i], label="dim", color=COLOR1)
+        ax.fill_between(
+            t[i], y1[i] - y1_err[i], y1[i] + y1_err[i], color=COLOR3, alpha=0.3
+        )
+        ind = int(np.round(0.2 * len(t[i])))
+        angle = (
+            360
+            / (2 * np.pi)
+            * np.atan(
+                (y1[i][ind + 1] - y1[i][ind])
+                / (np.max(y1) - np.min(y1))
+                / (t[i][ind + 1] - t[i][ind])
+                * (np.max(t) - np.min(t))
+            )
+        )
+        y = y1[i][ind] + 0.1 * (np.max(y1[i]) - np.min(y1[i]))
+        ax.text(t[i][ind], y, f"D={diffusion_constants[i]}", rotation=angle)
     ax.set_ylabel("Fractal Dimension")
     ax.set_xlabel("Time [min]")