Analyzed scaling of cell lists

gen/cell_scaling.py

+import sys
+sys.path.append('../')
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from netflow import *
+
+def main():
+    # Load basnet
+    basenet = netflow.load_network_from('../netflow/network/network.h5')
+    print("Network statistics:")
+    print("Throats: %d" % len(basenet.throats))
+    pore_throat_counts = [len(pore.throats) for pore in basenet.pores]
+    print("Max. number of throats per pore: %d" % max(pore_throat_counts))
+    print("Min. number of throats per pore: %d" % min(pore_throat_counts))
+    print("Avg. number of throats per pore: %f" % (sum(pore_throat_counts) / len(pore_throat_counts)))
+    
+    cutoff  = 0.5 * max([basenet.ub[i] - basenet.lb[i] \
+                        for i in range(len(basenet.ub))])
+    
+    end = 8
+    target_sizes = np.arange(end) + 1
+    avg_vals = np.empty(end, dtype=np.float64)
+    max_vals = np.empty(end, dtype=np.float64)
+    for i in range(1, end+1):
+        target = [i, i, i]
+        print(target)
+        max_, _, avg = netflow.cell_list_scaling(basenet=basenet, targetsize=target, \
+                                  cutoff=cutoff, sd=42)
+        avg_vals[i-1] = avg
+        max_vals[i-1] = max_
+    # Plot scaling of average number of pores per cell for different targetsizes
+    plt.figure()
+    plt.title(r"Scaling of Cell-List statistics w.r.t. Target-Size")
+    plt.xlabel(r"target size")
+    plt.ylabel(r"average #pores per cell")
+    plt.plot(target_sizes, avg_vals, "--bo")
+    plt.grid()
+    plt.savefig("plots/cell_scaling_avg.png")
+    plt.close()
+
+    plt.figure()
+    plt.title(r"Scaling of Cell-List statistics w.r.t. Target-Size")
+    plt.xlabel(r"target size")
+    plt.ylabel(r"maximum #pores per cell")
+    plt.plot(target_sizes, max_vals, "--ro")
+    plt.grid()
+    plt.savefig("plots/cell_scaling_max.png")
+    plt.close()
+    
+if __name__ == '__main__':
+    main()

gen/generate.py

@@ -10,9 +10,10 @@ def main():
     basenet = netflow.load_network_from('../netflow/network/network.h5')
     print("Network statistics:")
     print("Throats: %d" % len(basenet.throats))
-    print("Max. number of throats per pore: %d" % max([len(pore.throats) for pore in basenet.pores]))
-    print("Min. number of throats per pore: %d" % min([len(pore.throats) for pore in basenet.pores]))
-    print("Avg. number of throats per pore: %f" % (sum([len(pore.throats) for pore in basenet.pores]) / len(basenet.pores)))
+    pore_throat_counts = [len(pore.throats) for pore in basenet.pores]
+    print("Max. number of throats per pore: %d" % max(pore_throat_counts))
+    print("Min. number of throats per pore: %d" % min(pore_throat_counts))
+    print("Avg. number of throats per pore: %f" % (sum(pore_throat_counts) / len(pore_throat_counts)))
     
     target  = [3, 3, 3]
     cutoff  = 0.5 * max([basenet.ub[i] - basenet.lb[i] \

netflow/netgen.py

@@ -120,22 +120,6 @@ class CellList:
         for pore in pores:
             cellIdx = self.pore_to_index(pore)
             self.poresSorted[cellIdx].add(pore)
-        """ DEBUG
-        cellCounts = [0] * self.totalCells
-        for pore in pores:
-            cellIdx = self.pore_to_index(pore)
-            cellCounts[cellIdx] += 1
-        n = len(pores)
-        max_ = max(cellCounts)
-        min_ = min(cellCounts)
-        avg  = sum(cellCounts) / len(cellCounts)
-        print("Cell statistics:")
-        print("Number of cells (%d, %d, %d)" % (self.nCells[0], self.nCells[1], self.nCells[2]))
-        print("Max. pores per cell: %d percent: %f" % (max_, max_ / n * 100.))
-        print("Min. pores per cell: %d percent: %f" % (min_, min_ / n * 100.))
-        print("Avg. pores per cell: %f percent: %f" % (avg, avg / n * 100.))
-        raise AssertionError("Stop")
-        """
         
     # Helper for converting 3D index triple into 1D (flattened) index
     def flatten(self, i: int, j: int, k: int) -> int:
@@ -180,6 +164,22 @@ class CellList:
                             candidates[nbor] = l
         return candidates
     
+    def cell_stats(self):
+        cellCounts = [len(self.poresSorted[i]) for i in range(self.totalCells)]
+        n    = sum(cellCounts)
+        max_ = max(cellCounts)
+        min_ = min(cellCounts)
+        avg  = n / self.totalCells
+        
+        """
+        print("Cell statistics:")
+        print("Number of cells (%d, %d, %d)" % (self.nCells[0], self.nCells[1], self.nCells[2]))
+        print("Max. pores per cell: %d (%f%%)" % (max_, max_ / n * 100.))
+        print("Min. pores per cell: %d (%f%%)" % (min_, min_ / n * 100.))
+        print("Avg. pores per cell: %.1f (%f%%)" % (avg, avg / n * 100.))
+        """
+        return (max_, min_, avg)
+
     # Remove pore from sorted pores
     def expel(self, pore: Pore):
         poreIdx = self.pore_to_index(pore)
@@ -454,6 +454,98 @@ def generate_imperial(basenet: Network, targetsize: List[float], \
     for throat in throats: network.connect_pores(pore1=pores[throat[0]],
         pore2=pores[throat[1]], label=throat[2], r=throat[3])
     return network
+
+def cell_list_scaling(basenet: Network, targetsize: List[int], \
+    cutoff: float = float('inf'), sd: int = None, mute: bool = False):
+    from random import seed, random, randint
+    from itertools import product
+    seed(sd)
+    d = len(targetsize) # number of spatial dimensions
+    # check target size multiplicator (must be >= 1)
+    if any([i < 1 for i in targetsize]):
+        raise NameError('targetsize must be >= 1!')
+    # size of new network
+    L = list(map(lambda lb,ub,i: (ub-lb)*i, basenet.lb, basenet.ub, targetsize))
+    # check target size (must be > basenet.Lmax)
+    if any([Lc < basenet.Lmax for Lc in L]):
+        raise NameError('targetsize leads network < Lmax of basenet!')
+
+    # pores, distributed based on dendrogram of basenet
+    if (not mute): print("distributing pores...", end="", flush=True)
+    # make indexable and discard in-/outflow pores
+    basepores = [pore for pore in basenet.pores \
+        if ((pore.label[:len(LABELS[1])] != LABELS[1]) \
+        and (pore.label[:len(LABELS[2])] != LABELS[2]))]
+    # dendrogram-based or uniform pore distribution
+    pores = []
+    if (cutoff != cutoff): # uniform
+        print("\b"*21 + "uniform pore distribution")
+        # loop over network sections & add uniformly distributed pores drawn from basenet
+        for s in product(*[range(k) for k in targetsize]):
+            for k, pore in enumerate(basepores):
+                pos = [random()*l for l in L]
+                pores.append(Pore(pos=pos, r=pore.r, label=LABELS[0], 
+                    throats=pore.throats.copy()))
+    else: # dendrogram-based
+        # extract cluster hierarchy from basenet
+        centroids = [pore.pos for pore in basepores]
+        from scipy.cluster import hierarchy
+        clustree = hierarchy.linkage(centroids, method = 'centroid')
+        # determine cluster centroids and weights
+        weights = len(basepores)*[1] # pores have weight = 1
+        for cluster in clustree:
+            p1 = int(cluster[0]); p2 = int(cluster[1])
+            w1 = weights[p1]; w2 = weights[p2]
+            centroids.append(list(map(lambda x1,x2: (w1*x1 + w2*x2)/(w1 + w2),
+                centroids[p1], centroids[p2])))
+            weights.append(w1 + w2)
+        # loop over network sections (twisted basenet copies)
+        for s in product(*[range(k) for k in targetsize]):
+            # twist centers of sufficiently small cluster
+            touched = [False]*len(centroids)
+            # positions of pores based on rotations of linked cluster-pairs
+            for k in range(len(centroids)-1,len(basepores)-1,-1):
+                i = k-len(basepores) # index of cluster in clustree
+                ctr = centroids[k] # center of rotation
+                dist = clustree[i,2] # cluster distance
+                # check if cluster/pore needs to be twisted
+                if ((dist < cutoff) or touched[k]):
+                    p1 = int(clustree[i,0]); p2 = int(clustree[i,1])
+                    w1 = weights[p1]; w2 = weights[p2]
+                    vec = random_direction(d)
+                    centroids[p1] = \
+                        [ctr[j] + w2*dist/(w1+w2)*vec[j] for j in range(d)]
+                    centroids[p2] = \
+                        [ctr[j] - w1*dist/(w1+w2)*vec[j] for j in range(d)]
+                    touched[p1] = touched[p2] = True
+            # add section pores in random order to pores of new network
+            for k, pore in enumerate(basepores):
+                pos = [c-lb + si*(ub-lb) \
+                    for c,si,lb,ub in zip(centroids[k],s,basenet.lb,basenet.ub)]
+                pores.insert(randint(0,len(pores)), Pore(pos=pos,
+                    r=pore.r, label=LABELS[0], throats=pore.throats.copy()))
+        print("\b"*21 + "left {0:d} of {1:d} clusters (incl. {2:d} pores) untouched".\
+            format(sum([int(not j) for j in touched]), len(touched), len(basepores)))
+        # flip pores outside back into domain
+        for pore in pores:
+            for k in range(d):
+                if pore.pos[k] < 0:
+                    pore.pos[k] = L[k] + pore.pos[k]
+                elif pore.pos[k] >= L[k]:
+                    pore.pos[k] = pore.pos[k] - L[k]
+
+    # add pore buffer layers for spatial periodicity
+    n = len(pores)
+    copies = __add_buffer_layers(pores, L, basenet.Lmax)
+    
+    # Domain size including periodic buffer layers on all sides
+    trueDomainSize = [L[i] + 2 * basenet.Lmax for i in range(d)]
+    # Initialize cell-list; Place each pore in respective cell-set
+    cellList = CellList(pores, trueDomainSize, basenet.Lmax, basenet.Lmax)
+    # Print cell stats to console
+    max_, min_, avg = cellList.cell_stats()
+    return (max_, min_, avg)
+    
     
 def generate_dendrogram(basenet: Network, targetsize: List[int], \
     cutoff: float = float('inf'), sd: int = None, mute: bool = False) \