netgen.py

# generation and planar cutting of periodic pore networks
#
# Daniel W. Meyer
# Institute of Fluid Dynamics, ETH Zurich
# January 2019
import multiprocessing as mp
from typing import List, Dict, Set, Tuple # for type hints in argument lists

LABELS = ('', 'in', 'out', 'cut') # don't change the order

class Pore:
    """Class of a pore connected to other pores via throats."""
    id = 0
    def __init__(self, pos: List[float], r: float, label: str = LABELS[0], \
        throats: Set = None, id: int = -1, index: int = -1):
        if id == -1: # no id given
            self.id = Pore.id; Pore.id = Pore.id+1
        else: # id provided
            self.id = id; Pore.id = max(id+1, Pore.id)
        self.pos = pos.copy() # position vector
        if throats is None: throats = set()
        self.throats = throats # throats set
        self.r = r # radius
        self.label = label # label string like '', 'in', or 'out'
        # Identifier for pores in list used in generate_dendrogram
        self.index = index
    def __repr__(self): return str(self.__class__) + ': ' + \
        str({k:self.__dict__[k] for k in self.__dict__ if k != 'throats'}) + \
        ' {0:d} throats'.format(len(self.throats)) # dont flush throat objects

class Throat:
    """Class of a throat connecting two pores.
    
    Periodic throats have a label 'X1 X2 X3' with Xc being element of {-1,0,1}.
    For Xc = 1, pore1 is at the right or upper domain bound in the c-direction
    and pore2 is at the left or lower bound. Vice versa for Xc = -1. With
    Xc = 0, the throat does not cross the bound periodically in c-direction.
    """
    id = 0
    def __init__(self, pore1: Pore, pore2: Pore, r: float, \
        label: str = LABELS[0], id: int = -1, index: int = -1):
        if id == -1: # no id given
            self.id = Throat.id; Throat.id = Throat.id+1
        else: # id provided
            self.id = id; Throat.id = max(id+1, Throat.id)
        # pore objects
        self.pore1 = pore1
        self.pore2 = pore2
        # radius and label
        self.r = r
        self.label = label # label string like '' or periodicity label '1 0 -1'
    def __repr__(self): return str(self.__class__) + ': ' + str(self.__dict__)

class Network:
    """Network class with pore and throat lists.
    
    Lmax must be smaller than the smallest side length of the network ub-lb.
    """
    def __init__(self, lb: List[float] = [], ub: List[float] = [], \
        pores: List[Pore] = None, throats: List[Throat] = None, \
        Lmax: float = 0.0, label: str = None):
        # lower/upper bounds for network volume
        self.lb = lb.copy(); self.ub = ub.copy()
        # sets with pores and throats
        if pores is None: pores = set()
        if throats is None: throats = set()
        self.pores = pores.copy(); self.throats = throats.copy()
        # length of longest throat (distance between connected pores)
        self.Lmax = Lmax
        # network label or name
        from datetime import datetime
        if label is None:
            self.label = datetime.today().isoformat(sep=' ',timespec='seconds')
        else:
            self.label = label
    def __repr__(self):
        # count in/out resp. labelled pores
        k = sum([int(pore.label != LABELS[0]) for pore in self.pores])
        # report
        return 'Network \'' + self.label + '\' from ' + str(self.lb) + \
            ' to ' + str(self.ub) + \
            ' with {0:d}({1:d}) pores, {2:d} throats, Lmax = {3:e}'. \
            format(len(self.pores), k, len(self.throats), self.Lmax)
    def add_pore(self, pore: Pore):
        self.pores.add(pore)
    def remove_pore(self, pore: Pore):
        """Remove pore and all throats connected to it."""
        # disconnect throats from pores connected to pore
        for throat in pore.throats.copy():
            connected_pore = throat.pore1
            if (connected_pore == pore): connected_pore = throat.pore2
            connected_pore.throats.remove(throat)
            pore.throats.remove(throat)
            self.throats.remove(throat)
        # remove pore
        self.pores.remove(pore)
    def connect_pores(self, pore1: Pore, pore2: Pore, r: float, \
        throat_id: int = -1, label: str = LABELS[0]) -> Throat:
        """Establish a throat connection between two pores."""
        throat = Throat(pore1, pore2, r, label, throat_id)
        self.throats.add(throat)
        pore1.throats.add(throat); pore2.throats.add(throat)
        return throat

class CellList:
    """CellList class dividing physical domain into equally-sized cells
    
    find_candidates() assumes pore is located in interior of domain (non-periodic)
    """
    def __init__(self, pores: List[Pore], copies: Dict[Pore, Set[Pore]],  \
                    domainSize: List[float], cellSize: float, \
                        basenet: Network, nInterior: int, isPeriodicCheck: bool):
        from math import ceil
        self.dim = len(domainSize)
        self.Lmax = basenet.Lmax
        # Number of cells in each dimension
        self.nCells = [max(1, ceil(domainSize[i] / cellSize)) for i in range(self.dim)]
        # Inverse of cell sizes in domain
        self.invCellSizes = [self.nCells[i] / domainSize[i] for i in range(self.dim)]
        # Total number of cells
        self.totalCells = prod(self.nCells)
        # Pore positions 
        self.poresPos = []  # size: dim * #pores
        # Throat lengths
        self.throatL = []  # size: #throats (from non-periodic pores)
        # Throat radii
        self.throatR = []  # size: #throats (from non-periodic pores)
        throatIdx = 0
        # Offset in pore-throat table for given pore index
        self.throatOffsets = [0]  # size: #pores + 1 (non-periodic)
        offset = 0
        # Used to keep track of already realized throats
        self.poreThroatIndices = [set() for _ in range(nInterior)]
        # Additional information needed for correctness if targetsize contains 1
        self.isPeriodicCheck = isPeriodicCheck
        if isPeriodicCheck:
            # Maps periodic copy index to respective index inside domain
            # {int: int}
            self.originalIdx = {}
            # Maps index from inside domain to List of all its periodic copies
            # {int: List[int]}
            self.periodicCopies = {}
            
        L = [ub-lb for lb,ub in zip(basenet.lb,basenet.ub)]
        # Sort pores according to cell membership
        self.poresSorted = [set() for _ in range(self.totalCells)]
        
        if isPeriodicCheck:
            for poreIdx, pore in enumerate(pores):
                cellIdx = self.pore_to_index(pore)
                
                self.poresSorted[cellIdx].add(poreIdx)
                self.poresPos += pore.pos
                
                if poreIdx < nInterior:
                    # Pre-fix sum
                    nthroats = len(pore.throats)
                    
                    offset += nthroats
                    self.throatOffsets.append(offset)
                    
                    self.originalIdx[poreIdx] = poreIdx  # identity
                    
                    copyIndices = []
                    for cpore in copies[pore]:
                        copyIndices.append(cpore.index)
                    self.periodicCopies[poreIdx] = copyIndices
                        
                    for throat in pore.throats:
                        self.throatL.append(distance(*throat_ends(throat, L)))
                        self.throatR.append(throat.r)
                        self.poreThroatIndices[poreIdx].add(throatIdx)
                        
                        throatIdx += 1
                else:  # periodic copy of pore in interior
                    # Determine index of original pore
                    originalIdx, _ = pore.label.split(' ',1)
                    self.originalIdx[poreIdx] = int(originalIdx)
        
        else:
            for poreIdx, pore in enumerate(pores):
                cellIdx = self.pore_to_index(pore)
                
                self.poresSorted[cellIdx].add(poreIdx)
                self.poresPos += pore.pos
                
                if poreIdx < nInterior:
                    # Pre-fix sum
                    nthroats = len(pore.throats)
                    
                    offset += nthroats
                    self.throatOffsets.append(offset)
                        
                    for throat in pore.throats:
                        self.throatL.append(distance(*throat_ends(throat, L)))
                        self.throatR.append(throat.r)
                        self.poreThroatIndices[poreIdx].add(throatIdx)
                        
                        throatIdx += 1
        
        throatTotal = offset
        self.poreMatchTable = mp.RawArray('i', throatTotal)
        
    # Helper for converting 3D index triple into 1D (flattened) index
    def flatten(self, i: int, j: int, k: int) -> int:
        return i + self.nCells[0] * (j + self.nCells[1] * k)
    
    def pore_to_triplet(self, poreIdx: int) -> List[int]:
        from math import floor
        # Shift pore position to be positive first (buffer layer)
        porePos = self.fetch_pos(poreIdx)
        return [floor( (porePos[i] + self.Lmax) * self.invCellSizes[i]) \
                        for i in range(self.dim)]
    # Compute index of given pore in cell list based on position
    def pore_to_index(self, pore: Pore) -> int:
        from math import floor
        # Shift pore position to be positive first (buffer layer)
        cellIdx = self.flatten(*[floor( (pore.pos[i] + self.Lmax) * self.invCellSizes[i]) \
                        for i in range(self.dim)])
        return cellIdx
    
    def is_valid_index(self, idx: int) -> bool:
        return 0 <= idx < self.totalCells
    
    def match_maker(self, pores: List[int], tid: int):
        n = len(pores)
        
        for poreIdx in pores:
            
            if tid == 0:
                print(f"progress {((poreIdx+1) / n)*100.:.1f}%", end="\r", flush=True)
                
            candidates = self.find_candidates(poreIdx)
            
            throats = self.fetch_throat_indices(poreIdx)
            for idx, throatIdx in enumerate(throats):
                lt = self.throatL[throatIdx]
                
                # Find best match among candidates
                matchIdx, matchLen = min(candidates.items(), key=lambda il: abs(il[1] - lt))
                # Remove match from candidates
                del candidates[matchIdx]
                # Put best match index in poreThroatTable
                self.set_throat(poreIdx, matchIdx, idx)
                # Remove periodic copies as well if targetsize contains 1
                if self.isPeriodicCheck:
                    # Get index of original pore in interior of domain
                    originalIdx = self.originalIdx[matchIdx]
                    
                    if originalIdx in candidates:
                        del candidates[originalIdx]
                    
                    for cporeIdx in self.periodicCopies[originalIdx]:
                        if cporeIdx in candidates:
                            del candidates[cporeIdx]
    
    # Compute dictionary of all nbor candidates of given pore
    def find_candidates(self, poreIdx: int) -> Dict[int, float]:    
        candidates  = {}
        
        porePos = self.fetch_pos(poreIdx)
        neighborhood = self.nbor_indices(poreIdx)
        nearestCellIdx = neighborhood[0]
        
        # Search current cell
        for nborIdx in self.poresSorted[nearestCellIdx]:
            if (nborIdx == poreIdx): continue
            l = distance(self.fetch_pos(nborIdx), porePos)
            if (l < self.Lmax):
                candidates[nborIdx] = l
        # Search neighboring cells
        for nborCellIdx in neighborhood[1:]:
            for nborIdx in self.poresSorted[nborCellIdx]:
                l = distance(self.fetch_pos(nborIdx), porePos)
                if (l < self.Lmax):
                    candidates[nborIdx] = l
        
        return candidates
    
    def set_throat(self, poreIdx: int, matchIdx: int, idx: int):
        index = self.throatOffsets[poreIdx] + idx
        # Set matchIdx
        self.poreMatchTable[index] = matchIdx
        
    def fetch_matches(self, poreIdx: int) -> List[int]:
        lb = self.throatOffsets[poreIdx]
        ub = self.throatOffsets[poreIdx + 1]
        return self.poreMatchTable[lb:ub]
    
    def fetch_throat_indices(self, poreIdx: int) -> List[int]:
        lb = self.throatOffsets[poreIdx]
        ub = self.throatOffsets[poreIdx + 1]
        return range(lb, ub)
        
    def fetch_pos(self, poreIdx: int) -> List[float]:
        lb = self.dim * poreIdx
        ub = self.dim + lb
        return self.poresPos[lb:ub]
    
    # Only consider subset of all neighboring cells (closest 7)
    def nbor_indices(self, poreIdx: int) -> List[int]:
        cellIdx = self.pore_to_triplet(poreIdx)
        # i, j, k; i +- 1, j, k; i, j +- 1, k; i, j, k +- 1;
        base   = self.flatten(*cellIdx)
        widthx = self.nCells[0]
        widthy = widthx * self.nCells[1]
        return [base, base - 1, base + 1, \
                base - widthx, base + widthx, \
                base - widthy, base + widthy]
        
    def nbor_cell_counts(self, neighborhood: List[int]) -> List[int]:
        return [len(self.poresSorted[i]) for i in neighborhood]
        
    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):
        poreCellIdx = self.pore_to_index(pore)
        self.poresSorted[poreCellIdx].remove(pore.index)

def distance(p1: List[float], p2: List[float]):
    """Distance between points p1 and p2."""
    mag = 0.0
    for c1, c2 in zip(p1,p2): mag += (c1-c2)**2
    return mag**0.5


def prod(v):
    """Cumulative product of vector components."""
    from functools import reduce
    return(reduce(lambda a,b: a*b, v))


def random_direction(d: int) -> List[float]:
    """Unity vector with uniformly distributed orientation."""
    from math import pi, sin, cos, acos
    from random import random
    phi = 2*pi * random()
    if (d == 2):
        return([cos(phi), sin(phi)])
    elif (d == 3):
        theta = acos(2*random() - 1)
        return([sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)])
    else:
        raise ValueError("random_direction supports d = 2 or 3")


def throat_ends(throat: Throat, L: List[float]) -> Tuple[List[float]]:
    """Provide the correct relative end points of a possibly periodic throat."""
    x1 = throat.pore1.pos.copy(); x2 = throat.pore2.pos.copy()
    # account for periodic throats
    if len(throat.label) != 0:
        direct = [int(k) for k in throat.label.split()] # periodicity
        x2 = [x2[k] + direct[k] * L[k] for k in range(len(L))]
    return (x1, x2)



def imperial_read(net_file_pref: str) -> Network:
    """Read network data in Imperial College format.
    
    For format specifics see PhD thesis of Taha Sochi from 2007 at Imperial.
    The names of the network files starts with the string net_file_pref.
    A Network object is returned.
    """
    # initialize network
    from os import sep
    network = Network(label = net_file_pref.split(sep)[-1])
    
    # read pore data
    with open(net_file_pref + '_node1.dat', 'r') as file:
        line = file.readline()
        size = [float(n) for n in line.replace('\t',' ').split()][1:]
        d = len(size) # dimensionality
        f = lambda words: (
            int(words[0]), # id
            [float(n) for n in words[1:(d+1)]], # pos
            [int(throat) for throat in words[(d+4+int(words[d+1])):]], # throats
            LABELS[int(words[d+2+int(words[d+1])]) + \
            2*int(words[d+3+int(words[d+1])])]) # label
        pores = [f(line.replace('\t',' ').split()) for line in file]
    pores.sort(key=lambda pore: pore[0])
    pores = [Pore(pore[1], 0.0, pore[3]) for pore in pores]
    # read pore radii
    with open(net_file_pref + '_node2.dat', 'r') as file:
        f = lambda words: (int(words[0]), float(words[2]))
        radii = [f(line.replace('\t',' ').split()) for line in file]
    radii.sort(key=lambda pore: pore[0])
    for k, pore in enumerate(pores): pore.r = radii[k][1]
    # add pores to network
    for pore in pores: network.add_pore(pore)

    # read throat data
    with open(net_file_pref + '_link1.dat', 'r') as file:
        n = int(file.readline().strip('\n ')) # number of throats
        f = lambda words: (
            int(words[0]), int(words[1])-1, int(words[2])-1, float(words[3]))
        throats = [f(line.replace('\t',' ').split()) for line in file]
    throats.sort(key=lambda throat: throat[0])
    # connect pores
    for throat in throats:
        if (throat[1] >= 0) and (throat[2] >= 0): # no in-/outflow throats
            network.connect_pores(pores[throat[1]], pores[throat[2]], throat[3])

    # determine maximum throat length
    Lmax = 0.0
    for throat in network.throats:
        L = distance(throat.pore1.pos, throat.pore2.pos)
        Lmax = max(L,Lmax)
    network.Lmax = Lmax

    network.lb = [0 for i in range(len(size))]
    network.ub = size
    return network


def plot_network(network: Network, labels: bool = False):
    """Plot pore network.
    
    Returns a figure object. If labels == True, pore and throat ids are shown.
    Pores with label == '' are plotted in black. Pores with labels starting
    with 'in', 'out', 'cut' are plotted in red, blue, magenta, respectively.
    """
    # setup plot
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    plt.ion()
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    
    # plot throats
    for throat in network.throats:
        if throat.label == LABELS[0]: col = 'black'
        else: col = 'grey'
        p1 = throat.pore1.pos; p2 = throat.pore2.pos
        ax.plot([p1[0], p2[0]], [p1[1], p2[1]], [p1[2], p2[2]],
            linewidth=0.5, color=col)
        if labels:
            ax.text3D((p1[0]+p2[0])/2, (p1[1]+p2[1])/2, (p1[2]+p2[2])/2,
                str(throat.id), fontsize=5)

    # plotting function for pores
    def plot_pores(label, size, color):
        pnts = [pore.pos for pore in network.pores \
            if (pore.label[:max(1,len(label))] == label)]
        if len(pnts) == 0: return # in case of empty point list
        pnts = list(map(list, zip(*pnts))) # transpose double list
        ax.scatter(*pnts, s=size, marker='.', c=color)
        if labels:
            for pore in network.pores:
                if (pore.label[:max(1,len(label))] == label):
                    ax.text3D(*pore.pos, str(pore.id), fontsize=7, color=color)
    # plot ordinary pores
    plot_pores(LABELS[0], 4, 'black')
    # plot inflow pores
    plot_pores(LABELS[1], 7, 'red')
    # plot outflow pores
    plot_pores(LABELS[2], 7, 'blue')
    # plot cut pores
    plot_pores(LABELS[3], 7, 'magenta')

    #ax.axis('equal')
    return fig


def generate_simple_net(n_pores: int, targetsize: List[float], \
    r_pore: float, r_throat: float, coordinatnumb: int, Lmax: float, \
    sd: int = None) -> Network:
    """Generate a spatially periodic network with uniform pore distribution."""
    d = len(targetsize) # number of spatial dimensions
    density = n_pores / prod(targetsize)
    origin = [0.0 for k in range(d)]
    Lb = (1/density)**(1/d)
    basenet = Network(origin, [Lb for k in range(d)], \
        [Pore(origin, r_pore, throats = {k for k in range(coordinatnumb)})], \
        [Throat(None, None, r_throat)], Lmax)
    basenet.label = 'simplenet_' + basenet.label
    return generate_imperial(basenet, targetsize, sd, False)


def generate_imperial(basenet: Network, targetsize: List[float], \
    sd: int = None, correlated: bool = True) -> Network:
    """Generate a new spatially periodic network.

    Based on an existing network basenet, generate a new network of size
    targetsize by using the algorithm outlined on p.54 of the PhD thesis of
    Nasiru Abiodun Idowu from 2009 at Imperial College.
    Pores are uniformly distribute with pore number density as in basenet.
    Pore radii and number of throat connections are sampled from basenet.
    Closest pores are connected. Throat radii are taken from basenet such
    that throat radius and radius-sum of the connected pores are correlated.
    Periodic throats are marked with a periodicity label, e.g., '1 0 -1' in 3d.
    """
    from random import seed, random, randint
    seed(sd)
    d = len(targetsize) # number of spatial dimensions
    # check target size (must be > basenet.Lmax)
    if any([L < basenet.Lmax for L in targetsize]):
        raise NameError('targetsize must be > Lmax of basenet!')

    # pores, uniformly distributed
    # number of pores in new network
    import math
    basesize = list(map(lambda a,b: a-b, basenet.ub, basenet.lb))
    n = math.ceil(prod(targetsize) * len(basenet.pores) / prod(basesize))
    # distribute pores uniformly
    basepores = [pore for pore in basenet.pores] # make indexable
    pores = [Pore(pos=[random()*L for L in targetsize], # pos
        # radius, label, number of throat connections
        r=basepores[idx].r, label=LABELS[0],
        throats=len(basepores[idx].throats)) for j, idx in \
        enumerate([randint(0,len(basepores)-1) for k in range(n)])]

    # add pore buffer layers for spatial periodicity
    n = len(pores)
    _ = __add_buffer_layers(pores, targetsize, basenet.Lmax)

    # throats, connect pores
    # triangulation
    import numpy as np
    from scipy.spatial import Delaunay
    points = np.array([pore.pos for pore in pores])
    tri = Delaunay(points) # possibility for missing (coplanar) points!
    # return neighbors of point pnt based on triangulation tri
    def neighbors_of(pnt, tri):
        return tri.vertex_neighbor_vertices[1][ \
            tri.vertex_neighbor_vertices[0][pnt]: \
            tri.vertex_neighbor_vertices[0][pnt+1]]
    # connect pores
    throats = []
    for k, pore in enumerate(pores[:n]):
        # initialize neighbor lists
        nbors = neighbors_of(k, tri)
        nbors = dict(zip(nbors,
            [distance(pore.pos, pores[j].pos) for j in nbors]))
        oldies = {k:0.0}
        # establish throat connections between pore and its neighbors
        while (pore.throats > 0) \
            and (len(nbors) > 0) and (min(nbors.values()) <= basenet.Lmax):
            # find closest neighbor
            for nbor in nbors:
                if nbors[nbor] == min(nbors.values()): break
            # if nbor is periodic copy, get original pore and periodicity label
            if (pores[nbor].label != LABELS[0]):
                nbor_o, lbl = pores[nbor].label.split(' ',1)
                nbor_o = int(nbor_o)
            else:
                nbor_o = nbor; lbl = LABELS[0]
            # is nbor a valid pore to connect to?
            if not ((nbor_o <= k) or (nbor in oldies) or \
                (pores[nbor_o].throats == 0)):
                # yes -> add connection to pore nbor_o
                throats.append([k, nbor_o, lbl])
                pore.throats = pore.throats - 1
                pores[nbor_o].throats = pores[nbor_o].throats - 1
            # nbor has been dealt with
            oldies[nbor] = nbors[nbor]
            # expand search neighborhood
            for j in neighbors_of(nbor, tri):
                nbors[j] = distance(pores[k].pos, pores[j].pos)
            for oldie in oldies.keys():
                if oldie in nbors: del nbors[oldie]
    
    # initialize throat sets of pores
    k = 0 # count of unrealised throats
    for pore in pores[:n]:
        k = k + pore.throats
        pore.throats = set() # initialize throats set

    # set throat radii
    # sort branch weights (branch weight is prop. to sum of pore radii)
    bw = [pores[throat[0]].r + pores[throat[1]].r for throat in throats]
    bw = list(enumerate(bw))
    bw.sort(key=lambda w: w[1])
    # sample throat radii from basenet
    basethroats = [throat for throat in basenet.throats] # make indexable
    r = [randint(0,len(basethroats)-1) for k in range(len(bw))]
    r = [basethroats[k].r for k in r]
    if correlated: r = sorted(r)
    # assign radii
    for k in range(len(r)): throats[bw[k][0]].append(r[k])

    # assemble and return network
    network = Network(lb=[0.0 for k in range(d)],
        ub=targetsize, Lmax=basenet.Lmax, label='from_' + basenet.label)
    for pore in pores[:n]: network.add_pore(pore)
    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, nthreads: int = mp.cpu_count(), \
    mute: bool = False) -> Network:
    """Generate a new spatially periodic network.

    Based on an existing network basenet, generate a new network of size
    targetsize by using a dendrogram-based algorithm that accounts for the
    spatial distribution/clustering of pores.
    The pore-clustering dendrogram is taken from basenet. Pore radii and number
    of throats are sampled from basenet. Pores are connected based on the
    throat lengths from the basenet data. Throat radii are taken from basenet.
    Periodic throats are marked with a periodicity label, e.g., '1 0 -1' in 3d.
    If cutoff == float('nan'), pores are uniformly distributed.
    """
    import multiprocessing as mp
    from random import seed, random, randint
    from itertools import product, accumulate
    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...")
    # 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("uniform pore distribution")
        # loop over network sections & add uniformly distributed pores drawn from basenet
        runningIndex = 0
        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(), index=runningIndex))
                runningIndex += 1
    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("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 and set respective index of all pores
        for i, pore in enumerate(pores):
            pore.index = i
            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)
    
    # Flag for handling the case where original nbor pore and its periodic
    # copies are potential candidates of a given pore at the same time 
    # (necessary only if target size is 1 in at least one dimension)
    isPeriodicCheck = 1 in targetsize
    # 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, copies, trueDomainSize, basenet.Lmax, \
                            basenet, n, isPeriodicCheck)
    
    # Evenly distribute pores among threads
    remainder = n % nthreads
    loads = [n // nthreads] * nthreads
    for i in range(remainder):
        loads[i] += 1
    
    # Compute inclusive-scan
    displs = [0]
    displs += accumulate(loads)
    
    workers = []
    if (not mute): print("computing best matches")
    # Compute best matches for each pore in parallel using available threads
    for tid in range(nthreads):
        poreList = range(displs[tid], displs[tid+1])
        worker = mp.Process(target=cellList.match_maker, args=(poreList, tid))
        worker.start()
        workers.append(worker)
    
    for worker in workers:
        worker.join()
        
    # throats, connect pores
    if (not mute): print("\nconnecting")
    throats = []
    k = 0 # count of unrealised throats
    
    avg_throat_diff = 0.
    count = 0
    total = 0
    n_already_taken = 0
    for poreIdx, pore in enumerate(pores[:n]):
        # Delete throats (not needed anymore)
        pore.throats.clear()
        
        if (not mute):
            print(f"progress {((poreIdx+1) / n)*100.:.1f}%", end="\r", flush=True)
        # Set of throats to be realized
        ttbr = cellList.poreThroatIndices[poreIdx]
        if (len(ttbr) == 0): continue  # Nothing left to do
        
        matches = cellList.fetch_matches(poreIdx)
        throatIndices = cellList.fetch_throat_indices(poreIdx)
        
        for i, match in enumerate(matches):
            throatIdx = throatIndices[i]
            
            if (throatIdx not in ttbr):
                continue  # Throat already realized
                            
            nbor = pores[match]
            
            nborc = nbor # memorize potential periodic copy
            # find original of buffer layer pore
            if (nbor.label != LABELS[0]):
                j, lbl = nbor.label.split(' ',1)
                nbor = pores[int(j)]
            else:
                lbl = LABELS[0]
            
            nborTtbr = cellList.poreThroatIndices[nbor.index]
            
            total += 1
            if (len(nborTtbr) == 0):
                n_already_taken += 1
                continue
            
            r = cellList.throatR[throatIdx]
            targetLen = cellList.throatL[throatIdx]
            # Compute actual throat length between pores
            lt = distance(pore.pos, nborc.pos)
            
            # connect pore to nbor
            throats.append([pore, nbor, lbl, r])
            # remove most similar throat of nbor
            ranking = [(nborThroatIdx, abs(cellList.throatL[nborThroatIdx] - lt)) \
                         for nborThroatIdx in nborTtbr]
            nborThroatIdx, _ = min(ranking, key=lambda il: il[1])
            
            diff = abs(targetLen - lt)
            avg_throat_diff += diff
            count += 1
            
            # Remove throat from neighbor
            nborTtbr.remove(nborThroatIdx)
            # Successfully connected this throat
            ttbr.remove(throatIdx)
            
            # remove fully-connected nbor from search neighborhood
            if (len(nborTtbr) == 0):
                cellList.expel(nbor)
                
                for cpore in copies[nbor]:
                    cellList.expel(cpore)
        
        # remove fully-connected pore from search neighborhood
        if (len(ttbr) == 0):
            cellList.expel(pore)
        
            for cpore in copies[pore]:
                cellList.expel(cpore)
    
    # Free space
    del cellList.poreMatchTable
    
    # Finally, find next best candidates for left-out pores
    # in REVERSED order
    if (not mute): print("\nfinding next best candidates")
    
    poreIdx = n
    for it, pore in enumerate(reversed(pores[:n])):
        poreIdx -= 1
        
        if (not mute):
            print(f"progress {((it+1) / n)*100.:.1f}%", end="\r", flush=True)
        
        ttbr = cellList.poreThroatIndices[poreIdx]
        if (len(ttbr) == 0):
            continue  # Nothing left to do
        
        # Expensive: Find alternative (multiprocessing)
        candidates = cellList.find_candidates(poreIdx)
            
        for throatIdx in ttbr.copy():
            
            r = cellList.throatR[throatIdx]
            targetLen = cellList.throatL[throatIdx]
            
            if (len(candidates) == 0):
                # no candidates left; give up
                k += len(ttbr)  # update num unrealized throats
                break
                
            nborIdx, lt = min(candidates.items(), key=lambda il: abs(il[1] - targetLen))
                        
            nbor = pores[nborIdx]
            nborc = nbor # memorize potential periodic copy
            # find original of buffer layer pore
            if (nbor.label != LABELS[0]):
                j, lbl = nbor.label.split(' ',1)
                nbor = pores[int(j)]
            else:
                lbl = LABELS[0]
            
            # Remove nbor from candidates as well as all periodic copies
            del candidates[nborc.index]
            # Remove original non-periodic nbor as well as all its
            # periodic copies. Only applies if targetsize contains 1
            if isPeriodicCheck:
                if nbor.index in candidates:
                    del candidates[nbor.index]
                    
                for cpore in copies[nbor]:
                    if cpore.index in candidates:
                        del candidates[cpore.index]

            # connect pore to nbor
            throats.append([pore, nbor, lbl, r])
            
            nborTtbr = cellList.poreThroatIndices[nbor.index]
            # remove most similar throat of nbor
            ranking = [(nborThroatIdx, abs(cellList.throatL[nborThroatIdx] - lt)) \
                         for nborThroatIdx in nborTtbr]
            nborThroatIdx, _ = min(ranking, key=lambda il: il[1])
            
            diff = abs(targetLen - lt)
            avg_throat_diff += diff
            count += 1
            
            # Remove throat from pore
            ttbr.remove(throatIdx)
            # Remove throat from neighbor
            nborTtbr.remove(nborThroatIdx)
            
            # remove fully connected nbor from search neighborhood
            if (len(nborTtbr) == 0):
                cellList.expel(nbor)
                
                for cpore in copies[nbor]:
                    cellList.expel(cpore)
        
        # remove fully connected pore from search neighborhood
        if (len(ttbr) == 0):
            cellList.expel(pore)
        
            for cpore in copies[pore]:
                cellList.expel(cpore)
    
    # Free memory associated with cellList (not needed anymore)
    del cellList
    
    percent = k / (k + len(throats)) * 100.
    print("\b"*24 + "{0:d} throats in total, {1:d} unrealised ({2:.1f}%)".\
        format(len(throats), k, percent))