Pre assembly

lib/boxcount.py

@@ -54,6 +54,7 @@ df = read_xy_alt(
     "/Users/joshuacoles/Developer/checkouts/jc3091/CompB DLA/data-analysis/data/alpha/1.csv"
 )
 
+max_actual_radius = max(df.x.abs().max(), df.y.abs().max())
 max_radius = int(np.exp2(np.floor(np.log2(max_actual_radius)) + 1))
 data = np.zeros((max_radius * 2, max_radius * 2))
 data[df.x + max_radius, df.y + max_radius] = 1

lib/fd.py

@@ -0,0 +1,78 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Sep 16 09:47:15 2019
+@author: daniel
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def fractal_dimension(array, max_box_size=None, min_box_size=1, n_samples=20, n_offsets=0, plot=False):
+    """Calculates the fractal dimension of a 3D numpy array.
+
+    Args:
+        array (np.ndarray): The array to calculate the fractal dimension of.
+        max_box_size (int): The largest box size, given as the power of 2 so that
+                            2**max_box_size gives the sidelength of the largest box.
+        min_box_size (int): The smallest box size, given as the power of 2 so that
+                            2**min_box_size gives the sidelength of the smallest box.
+                            Default value 1.
+        n_samples (int): number of scales to measure over.
+        n_offsets (int): number of offsets to search over to find the smallest set N(s) to
+                       cover  all voxels>0.
+        plot (bool): set to true to see the analytical plot of a calculation.
+
+
+    """
+    # determine the scales to measure on
+    if max_box_size == None:
+        # default max size is the largest power of 2 that fits in the smallest dimension of the array:
+        max_box_size = int(np.floor(np.log2(np.min(array.shape))))
+    scales = np.floor(np.logspace(max_box_size, min_box_size, num=n_samples, base=2))
+    scales = np.unique(scales)  # remove duplicates that could occur as a result of the floor
+
+    # get the locations of all non-zero pixels
+    locs = np.where(array > 0)
+    voxels = np.array([(x, y, z) for x, y, z in zip(*locs)])
+
+    # count the minimum amount of boxes touched
+    Ns = []
+    # loop over all scales
+    for scale in scales:
+        touched = []
+        if n_offsets == 0:
+            offsets = [0]
+        else:
+            offsets = np.linspace(0, scale, n_offsets)
+        # search over all offsets
+        for offset in offsets:
+            bin_edges = [np.arange(0, i, scale) for i in array.shape]
+            bin_edges = [np.hstack([0 - offset, x + offset]) for x in bin_edges]
+            H1, e = np.histogramdd(voxels, bins=bin_edges)
+            touched.append(np.sum(H1 > 0))
+        Ns.append(touched)
+    Ns = np.array(Ns)
+
+    # From all sets N found, keep the smallest one at each scale
+    Ns = Ns.min(axis=1)
+
+    # Only keep scales at which Ns changed
+    scales = np.array([np.min(scales[Ns == x]) for x in np.unique(Ns)])
+
+    Ns = np.unique(Ns)
+    Ns = Ns[Ns > 0]
+    scales = scales[:len(Ns)]
+    # perform fit
+    coeffs = np.polyfit(np.log(1 / scales), np.log(Ns), 1)
+
+    # make plot
+    if plot:
+        fig, ax = plt.subplots(figsize=(8, 6))
+        ax.scatter(np.log(1 / scales), np.log(np.unique(Ns)), c="teal", label="Measured ratios")
+        ax.set_ylabel("$\log N(\epsilon)$")
+        ax.set_xlabel("$\log 1/ \epsilon$")
+        fitted_y_vals = np.polyval(coeffs, np.log(1 / scales))
+        ax.plot(np.log(1 / scales), fitted_y_vals, "k--", label=f"Fit: {np.round(coeffs[0], 3)}X+{coeffs[1]}")
+        ax.legend();
+    return (coeffs[0])