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, originIndex: 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'
        self.index = index
        self.originIndex = originIndex  # points to original pore inside domain
    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], domainSize: List[float], cellSize: float, \
                        basenet: Network, throatL: List[float], throatR: List[float], \
                        nInterior: int):
        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)
        # Offsets in throat arrays
        self.throatOffsets = [0]
        # Sort pores according to cell membership
        self.poresSorted = [set() for _ in range(self.totalCells)]
        # For each pore keep track of nbor pores connected to it
        self.connPores = [set() for _ in range(nInterior)]
        # List of (target) throat lengths from base network
        self.throatL = throatL
        # List of throat radii from base network
        self.throatR = throatR
        
        throatTotal = 0
        for i, pore in enumerate(pores):
            cellIdx = self.pore_to_index(pore)
            self.poresSorted[cellIdx].add(pore)
            
            if i < nInterior:
                throatTotal += len(pore.throats)
                self.throatOffsets.append(throatTotal)
                
        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, pore: Pore) -> List[int]:
        from math import floor
        # Shift pore position to be positive first (buffer layer)
        return [floor( (pore.pos[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
    
    # Called in parallel. Find all best neighbor pore candidates and
    # store their indices in poreMatchTable
    def match_maker(self, pores: List[Pore], tid: int):
        n = len(pores)
        for i, pore in enumerate(pores):
            if (tid == 0):
                print(f"progress {((i+1) / n)*100.:.1f}%", end="\r", flush=True)
            
            candidates = self.find_candidates(pore)
            # Find best candidate for each throat of pore
            for idx, throatIdx in enumerate(pore.throats):
                lt = self.throatL[throatIdx]
                
                if (len(candidates) == 0):
                    self.set_throat(pore.index, -1, idx)
                    break  # no possible candidates left
                # Find best match among candidates
                matchIdx, _ = min(candidates.items(), key=lambda il: abs(il[1] - lt))
                # Remove match from candidates
                del candidates[matchIdx]
                # Put best match index in poreMatchTable
                self.set_throat(pore.index, matchIdx, idx)
    
    # Compute dictionary of all nbor candidates of given pore
    def find_candidates(self, pore: Pore) -> Dict[int, float]:    
        candidates  = {}
        
        porePos = pore.pos
        neighborhood = self.nbor_indices(pore)
        nearestCellIdx = neighborhood[0]
        
        # Set of pore indices already connected to pore
        pConn = self.connPores[pore.index]
        # Search current cell
        for nbor in self.poresSorted[nearestCellIdx]:
            nborIdx = nbor.index
            # No self connection; Check if already connected to this pore
            if (nborIdx == pore.index or nbor.originIndex in pConn): continue
            l = distance(nbor.pos, porePos)
            if (l < self.Lmax):
                candidates[nborIdx] = l
        # Search neighboring cells
        for nborCellIdx in neighborhood[1:]:
            for nbor in self.poresSorted[nborCellIdx]:
                nborIdx = nbor.index
                # Check if already connected to this pore
                if (nbor.originIndex in pConn): continue
                l = distance(nbor.pos, porePos)
                if (l < self.Lmax):
                    candidates[nborIdx] = l
        
        return candidates
    
    # Set appropriate field in poreMatchTable to index of best candidate
    def set_throat(self, poreIdx: int, matchIdx: int, idx: int):
        index = self.throatOffsets[poreIdx] + idx
        # Set matchIdx
        self.poreMatchTable[index] = matchIdx
    
    # Retrieve indices of best matches for given pore after they have
    # been computed
    def fetch_matches(self, pore: Pore) -> List[int]:
        lb = self.throatOffsets[pore.index]
        ub = lb + len(pore.throats)
        return self.poreMatchTable[lb:ub]
    
    # Consider all 27 neighboring cells
    def nbor_indices(self, pore: Pore) -> List[int]:
        base   = self.pore_to_index(pore)
        widthx = self.nCells[0]
        widthy = widthx * self.nCells[1]
        return [base, base - 1, base + 1,
                base - widthx, base + widthx,
                base - widthy, base + widthy,
                base + widthx - 1, base + widthx + 1,
                base - widthx - 1, base - widthx + 1,
                base + widthy - 1, base + widthy + 1,
                base - widthy - 1, base - widthy + 1,
                base + widthx + widthy, base + widthx + widthy - 1,
                base + widthx + widthy + 1, base + widthx - widthy,
                base + widthx - widthy - 1, base + widthx - widthy + 1,
                base - widthx + widthy, base - widthx + widthy - 1,
                base - widthx + widthy + 1, base - widthx - widthy,
                base - widthx - widthy - 1, base - widthx - widthy + 1]
    
    def update_connection(self, poreIdx: int, nborIdx: int):
        self.connPores[poreIdx].add(nborIdx)
        self.connPores[nborIdx].add(poreIdx)
        
    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)
    
    def is_contained(self, pore: Pore) -> bool:
        poreCellIdx = self.pore_to_index(pore)
        return pore in self.poresSorted[poreCellIdx]
        
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!')
    # check valid number of threads
    if (nthreads < 1):
        raise ValueError('Number of threads must be >= 1!')
    
    # Initialize throat lengths and radii of base network to be used
    # later
    throatIdxMap = {}  # {int (throat id): int (throat index)}
    throatBaseIdx = 0
    throatL = []
    throatR = []
    baseL = [ub-lb for lb,ub in zip(basenet.lb,basenet.ub)]
    for pore in basenet.pores:
        for throat in pore.throats:
            lt = distance(*throat_ends(throat, baseL))
            throatL.append(lt)
            throatR.append(throat.r)
            throatIdxMap[throat] = throatBaseIdx
            throatBaseIdx += 1
        
    del baseL  # no longer needed
    # 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
        runningIndex = 0
        print("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]
                poreThroatIdx = set(throatIdxMap[throat] for throat in pore.throats)
                pores.append(Pore(pos=pos, r=pore.r, label=LABELS[0], 
                    throats=poreThroatIdx, 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
            baseThroatIdx = 0
            for k, pore in enumerate(basepores):
                [throat in pore.throats]
                pos = [c-lb + si*(ub-lb) \
                    for c,si,lb,ub in zip(centroids[k],s,basenet.lb,basenet.ub)]
                poreThroatIdx = set(throatIdxMap[throat] for throat in pore.throats)
                pores.insert(randint(0,len(pores)), Pore(pos=pos,
                    r=pore.r, label=LABELS[0], throats=poreThroatIdx))
        
        # Throat index map no longer needed
        del throatIdxMap
        
        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
            pore.originIndex = i  # pore is already the original
            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, \
                          throatL, throatR, n)
    
    # Divide by two (since we are counting every throat twice)
    totalThroats = len(cellList.poreMatchTable) // 2
    
    # throats, connect pores
    if (not mute): print("\nconnecting")
    throats = []
    
    # DEBUG
    """
    for i, pore in enumerate(pores):
        if (pore.label != LABELS[0]):
            j, _ = pore.label.split(' ',1)
            assert(pore.originIndex == int(j))
        assert(pore.index == i)
    """
    # coordination number for each pore
    #coords = [len(pore.throats) for pore in pores[:n]]
    avg_throat_diff = 0.
    count = 0
    
    # Throats-To-Be-Realized
    # Copy throats of pores to keep track of the ones already realized
    ttbrs = [pore.throats.copy() for pore in pores[:n]]
    
    if (not mute): print("Total throats: %d" % totalThroats)
    
    throatsLeft = totalThroats
    throatsUnrealized = totalThroats
    
    # Threshold of remaining pores at which we compute the best
    # matches serially instead of using threads
    serialThresh = nthreads * 4
    # Pores in interior that are not fully-connected yet
    poresRemain = pores[:n]
    while (throatsUnrealized > 0):
        
        if (not mute):
            print("Throats left: %d Throats Unrealized: %d (%.2f%%)" % \
                 (throatsLeft, throatsUnrealized, throatsLeft / totalThroats * 100.))
        
        nRemain = len(poresRemain)
        if (nRemain > serialThresh):
            # Compute loads and displacements for all threads
            remainder = nRemain % nthreads
            loads = [nRemain // nthreads + (i < remainder) for i in range(nthreads)]
        
            # Compute inclusive-scan yielding offsets in poresRemain threads
            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 = poresRemain[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()
        else:  # serial
            if (not mute): print("\nSearching serially...")
            cellList.match_maker(poresRemain, 0)
        
        # Iterate over remaining (reduced) pores
        for it, pore in enumerate(poresRemain):
            
            if (not mute):
                print(f"progress {((it+1) / n)*100.:.1f}%", end="\r", flush=True)
            # Set of throats to be realized
            poreTtbr = ttbrs[pore.index]
            if (len(poreTtbr) == 0): continue  # Nothing left to do
            # Set of indices of connected pores
            pConn = cellList.connPores[pore.index]
            matches = cellList.fetch_matches(pore)
            
            for throatIdx, match in zip(pore.throats, matches):
                
                if throatIdx not in poreTtbr:
                    continue  # Throat already realized
                
                if match == -1:  # give up
                    poreTtbr.clear()
                    break
                    
                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)]
                    assert(nborc.originIndex == int(j))
                else:
                    lbl = LABELS[0]
                
                nborTtbr = ttbrs[nbor.index]
                
                # Check that nbor is not fully connected and that
                # this pore is not already connected to it
                if (len(nborTtbr) == 0 or nbor.index in pConn):
                    continue
                # Update sets of connected pores
                cellList.update_connection(pore.index, nbor.index)
                
                r = throatR[throatIdx]
                targetLen = 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(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
                
                # Successfully connected this throat
                poreTtbr.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(poreTtbr) == 0):
                cellList.expel(pore)
            
                for cpore in copies[pore]:
                    cellList.expel(cpore)
        
        # DEBUG
        """
        for pore in poresRemain:
            matches = cellList.fetch_matches(pore)
            for match in matches:
                if (match == -1): break  # No match found
                mPore = pores[pores[match].originIndex]  # original pore
                assert(mPore in poresRemain)
        """
        # Copy throats back to original pores and update #throats left
        # and compute remaining (reduced) pores
        throatCount = 0
        nextPores = []
        for pore in poresRemain:
            poreTtbr = ttbrs[pore.index]
            pore.throats = poreTtbr.copy()
            
            if len(poreTtbr) > 0:
                throatCount += len(poreTtbr)
                # Insert in random order for improved load balancing
                nextPores.insert(randint(0,len(nextPores)), pore)
            else:
                # Free memory associated with conn-pores of fully
                # connected pores
                cellList.connPores[pore.index].clear()
        
        poresRemain = nextPores
        throatsUnrealized = throatCount // 2  # Counted all throats twice
        throatsLeft = totalThroats - len(throats)
        
    
    # Unrealized throats
    nUnrealized = throatsLeft
    
    # Free memory (not needed anymore)
    del cellList
    del ttbrs
    
    # DEBUG
    percent = nUnrealized  / totalThroats * 100.
    print("\b"*24 + "{0:d} throats in total, {1:d} unrealised ({2:.3f}%)".\
        format(len(throats), nUnrealized, percent))
    avg_throat_diff /= count
    print("Avg. throat length difference: %e" % avg_throat_diff)
    print("Relative to Lmax: %.1f%%" % (avg_throat_diff / basenet.Lmax * 100.))

    # assemble and return network
    network = Network(lb=[0.0 for k in range(d)],
        ub=L, Lmax=basenet.Lmax, label='from_' + basenet.label)
    for pore in pores[:n]:
        pore.throats.clear()  # remove (potentially) remaining pores
        network.add_pore(pore)
    for throat in throats: network.connect_pores(pore1=throat[0],
        pore2=throat[1], label=throat[2], r=throat[3])
    
    # DEBUG
    """
    # Check for correctness of generated network
    checkUnrealized = 0
    for pore in network.pores:
        nthroats = len(pore.throats)
        # Make sure pore coordination number is respected
        coordNum = coords[pore.index]
        assert(nthroats <= coordNum)
        checkUnrealized += coordNum - nthroats
        
        nborInd = set()
        for throat in pore.throats:
            assert(pore.index == throat.pore1.index or pore.index == throat.pore2.index)
            nbor = throat.pore2 if (throat.pore1.index == pore.index) else throat.pore1
            if (nbor.label != LABELS[0]):
                nbor = pores[int(j)]
            nborInd.add(nbor.index)
        # Make sure pore isn't connected to the same nbor twice
        assert(len(nborInd) == nthroats)
    
    checkUnrealized //= 2  # counted every throat twice
    # Make sure missing throats are accounted for
    assert(checkUnrealized == nUnrealized)
    """
    return network


def __add_buffer_layers(pores: List[Pore], targetsize: List[float],
    Lbuffer: float) -> Dict[Pore,Set[Pore]]:
    """Add pore buffer layers for spatial periodicity.
    
    Returns a dict where each key is a pore and the corresponding value
    is a set containing the copies of that pore.
    """
    n = len(pores) # number of pores inside domain
    d = len(targetsize) # number of spatial dimensions
    # label = pore index + periodicity label (-1,0,1) in each dim.
    for k, pore in enumerate(pores): pore.label = str(k) + d * ' 0'
    # add pore buffer layers
    copies = {pore:set() for pore in pores} # {pore:{set of periodic copies}}
    runningIndex = n
    for k in range(d):
        # left bound
        pores_layer = [pore for pore in pores \
            if targetsize[k]-Lbuffer <= pore.pos[k] < targetsize[k]]
        for pore in pores_layer:
            pcopy = Pore(pore.pos, pore.r, throats=pore.throats, \
                            id=pore.id, index=runningIndex)
            pcopy.pos[k] = pcopy.pos[k]-targetsize[k]
            lbl = pore.label.split(); lbl[k + 1] = '-1'
            pcopy.label = ' '.join(lbl) # periodicity label
            originIndex = int(lbl[0])
            pcopy.originIndex = originIndex
            pores.append(pcopy)
            copies[pores[originIndex]].add(pcopy)
            runningIndex += 1
        # right bound
        pores_layer = [pore for pore in pores \
            if 0 <= pore.pos[k] < Lbuffer]
        for pore in pores_layer:
            pcopy = Pore(pore.pos, pore.r, throats=pore.throats, \
                            id=pore.id, index=runningIndex)
            pcopy.pos[k] = pcopy.pos[k]+targetsize[k]
            lbl = pore.label.split(); lbl[k + 1] = '1'
            pcopy.label = ' '.join(lbl) # periodicity label
            originIndex = int(lbl[0])
            pcopy.originIndex = originIndex
            pores.append(pcopy)
            copies[pores[originIndex]].add(pcopy)
            runningIndex += 1
    # reset labels of pores inside domain (without buffer layers)
    for pore in pores[:n]: pore.label = LABELS[0]
    return copies


def open_periodic_network(network: Network, c: int, inouts: Dict = {}):
    """Open periodic throats in the spatial direction c>0 of a network.
    
    To avoid identical in/out pores resulting from different periodic throats
    connecting to the same original pore, one can keep track of newly generated
    in/out pores through the inouts pore dict {orig1:{"0 0 1":copy1, ...}, ...}.
    """
    # loop over (periodic) throats
    for throat in network.throats.copy(): # shallow copy (iteration changes set)
        # is throat a periodic connection in direction c?
        if len(throat.label) == 0: continue
        direct = [int(c) for c in throat.label.split()]
        if direct[c-1] == 0: continue
        # open periodic throat connection
        __open_periodic_throat(network, throat, direct, c, inouts)


def __open_periodic_throat(network: Network,
    throat: Throat, direct: List[int], c: int, inouts: Set[Pore]) \
    -> Set[Throat]:
    """Open a periodic throat connection.
    
    direct is the encoded periodicity label of the throat and c>0 is the
    spatial direction used for the in-/outflow labeling. The inouts dict
    is used to avoid the formation of pores with identical positions.
    """
    d = len(network.lb) # number of spatial dimensions
    p1 = throat.pore1; p2 = throat.pore2; r = throat.r
    # identify existing in/out pore or construct new one
    # needed in case of many periodic throats connecting to the same pore
    def get_inout_pore(pore_o, direct, pos, label, inouts) -> Pore:
        dirkey = str([abs(d) for d in direct]) # periodicity direction key
        # find existing inout pore copy (or periodic original)
        if pore_o in inouts: # is there an inout copy of the original pore?
            if dirkey in inouts[pore_o]: # is there copy in right direct?
                return inouts[pore_o][dirkey] # yes, return inout copy
        # if no inout copy exists, create a new one
        pore = Pore(pos, pore_o.r, label)
        if pore_o in inouts: # add new direct
            inout = inouts[pore_o] # dict dirkey:inout-pore
            inout[dirkey] = pore
        else: inouts[pore_o] = {dirkey:pore} # add original to inout dict
        return pore
    # pore2
    pos = [0.0 for j in range(d)] # init new position object
    # copy and shift pore2
    for j in range(d):
        pos[j] = p2.pos[j] + direct[j] * (network.ub[j]-network.lb[j])
    p2c = get_inout_pore(p2, direct, pos,
        LABELS[1+round((direct[c-1]+1)/2)] + str(c), inouts)
    network.add_pore(p2c)
    t1 = network.connect_pores(p1, p2c, throat.r)
    # pore1
    pos = [0.0 for j in range(d)] # init new position object
    # copy and shift pore1
    for j in range(d):
        pos[j] = p1.pos[j] - direct[j] * (network.ub[j]-network.lb[j])
    p1c = get_inout_pore(p1, direct, pos,
        LABELS[1+round((-direct[c-1]+1)/2)] + str(c), inouts)
    network.add_pore(p1c)
    t2 = network.connect_pores(p2, p1c, throat.r)
    # remove periodic throat
    network.throats.remove(throat)
    p1.throats.remove(throat); p2.throats.remove(throat)
    return {t1, t2}


def cut_network(network: Network, x: float, c: int, label: str = '',
    inouts: Dict = {}):
    """Cut a (periodic) network at position x_c with c>0.

    New pores are labelled and introduced at throat intersection points
    with the cutting plane. Before cutting, periodic throat connections
    in the c-direction are opened by calling open_periodic_network,
    which makes use of the inouts pore dictionary.
    """
    d = len(network.lb) # number of spatial dimensions
    # open periodic pores in direction c
    open_periodic_network(network, c, inouts)
    # cut throats incl. cut periodic ones being periodic normal to direction c
    throats = network.throats.copy() # shallow copy (iteration changes set)
    while len(throats) > 0:
        throats_next = set() # cut possible periodic pores in a 2nd round
        # cut throats ...
        for throat in throats:
            p1 = throat.pore1; p2 = throat.pore2
            # ... intersecting with cutting plane at x_c
            if (p1.pos[c-1] - x) * (x - p2.pos[c-1]) > 0:
                # treat periodic pores
                if len(throat.label) != 0:
                    direct = \
                        [int(k) for k in throat.label.split()] # periodicity
                    for k, cn in enumerate(direct):
                        if cn != 0: break # find first +/-1 periodic component
                    throats_periodic = __open_periodic_throat( \
                        network, throat, direct, k+1, inouts)
                    throats_next = throats_next | throats_periodic
                    continue
                # cut position
                f = (x - p1.pos[c-1]) / (p2.pos[c-1] - p1.pos[c-1])
                pos = [(p2.pos[j]-p1.pos[j])*f + p1.pos[j] for j in range(d)]
                pos[c-1] = x # to avoid tiny rounding errors
                # introduce and connect pore
                pore = Pore(pos, r=0.0, label=LABELS[3] + str(c) + label)
                network.add_pore(pore)
                # throat p1 - pore (reconnect throat from p2 to p1)
                throat.pore2 = pore; pore.throats.add(throat)
                # throat pore - p2 (new throat object)
                p2.throats.remove(throat)
                network.connect_pores(pore, p2, throat.r)
        # cut newly opened periodic throats
        throats = throats_next


def erase_network(network: Network, x: float, c: int, direct: bool, label: str):
    """Remove pores and connected throats from a network.
    
    Pores with positions < or > x_c for direct = True or False, respectively,
    are removed (c>0). At the same time, pores at x_c are labelled and network
    bounds updated.
    """
    for pore in network.pores.copy():
        # remove pores (and connected throats)
        if ((pore.pos[c-1] < x) and direct) or \
            ((x < pore.pos[c-1]) and not direct): network.remove_pore(pore)
        # relabel pores on cutting plane x_c
        elif pore.pos[c-1] == x: pore.label = label
    # update network bounds
    if direct: network.lb[c-1] = x
    else: network.ub[c-1] = x


def save_network_to(filename: str, network: Network):
    """Save pore network to hdf5 file."""
    import h5py
    f = h5py.File(filename, 'x')
    # write network label
    label = network.label.encode('ascii','ignore')
    f.create_dataset('network_label', (1,), dtype='S'+str(len(label)),
        data=[label])
    # bounds, Lmax
    f.create_dataset('lb', data=network.lb)
    f.create_dataset('ub', data=network.ub)
    f.create_dataset('Lmax', data=[network.Lmax])

    # write pore data
    p_grp = f.create_group('pores')
    # id
    wrk = [pore.id for pore in network.pores]
    p_grp.create_dataset('id', data=wrk)
    # radius
    wrk = [pore.r for pore in network.pores]
    p_grp.create_dataset('r', data=wrk)
    # position
    for k in range(len(network.lb)):
        wrk = [pore.pos[k] for pore in network.pores]
        p_grp.create_dataset('pos/x' + str(k), data=wrk)
    # label
    wrk = [pore.id for pore in network.pores if len(pore.label) != 0]
    if len(wrk) > 0:
        p_grp.create_dataset('label/id', data=wrk)
        wrk = [pore.label.encode('ascii','ignore') \
            for pore in network.pores if len(pore.label) != 0]
        from functools import reduce
        length = reduce(lambda a,b: max(a,b), [len(lbl) for lbl in wrk])
        p_grp.create_dataset('label/strg', (len(wrk),), dtype='S'+str(length),
            data=wrk)
        
    # write throat data
    t_grp = f.create_group('throats')
    # id
    wrk = [throat.id for throat in network.throats]
    t_grp.create_dataset('id', data=wrk)
    # radius
    wrk = [throat.r for throat in network.throats]
    t_grp.create_dataset('r', data=wrk)
    # pores
    wrk = [throat.pore1.id for throat in network.throats]
    t_grp.create_dataset('pore1', data=wrk)
    wrk = [throat.pore2.id for throat in network.throats]
    t_grp.create_dataset('pore2', data=wrk)
    # label
    wrk = [throat.id for throat in network.throats \
        if len(throat.label) != 0]
    if len(wrk) > 0:
        t_grp.create_dataset('label/id', data=wrk)
        wrk = [throat.label.encode('ascii','ignore') \
            for throat in network.throats if len(throat.label) != 0]
        from functools import reduce
        length = reduce(lambda a,b: max(a,b), [len(lbl) for lbl in wrk])
        t_grp.create_dataset('label/strg', (len(wrk),), dtype='S'+str(length),
            data=wrk)


def load_network_from(filename: str) -> Network:
    """Load pore network from hdf5 file."""
    import h5py
    import numpy as np
    f = h5py.File(filename, 'r')
    # global network properties
    network = Network(label=f['network_label'][0].decode('utf-8'),
        lb=list(f['lb']), ub=list(f['ub']), Lmax=f['Lmax'][0])

    # pores
    pid = np.array(f['pores/id'])
    pores = {} # id:pore dictionary
    for id in pid:
        pore = Pore(pos=[0.0 for k in network.lb], r=0.0, id=id)
        pores[id] = pore
        network.add_pore(pore)
    # radius
    for id, r in zip(pid, np.array(f['pores/r'])): pores[id].r = r
    # position
    for k in range(len(network.lb)):
        for id, x in zip(pid, np.array(f['pores/pos/x' + str(k)])):
            pores[id].pos[k] = x
    # labels
    if 'label' in f['pores']: # are there any labels != ''?
        for id, strg in zip(f['pores/label/id'],f['pores/label/strg']):
            pores[id].label = strg.decode('utf-8')
        
    # throats
    tid = np.array(f['throats/id'])
    throats = {} # id:throat dictionary
    for id, pore1, pore2, r in zip(tid,
        [pores[id] for id in np.array(f['throats/pore1'])],
        [pores[id] for id in np.array(f['throats/pore2'])],
        np.array(f['throats/r'])):
        throats[id] = network.connect_pores(pore1, pore2, r, id)

    # labels
    if 'label' in f['throats']: # are there any labels != ''?
        for id, strg in zip(f['throats/label/id'],f['throats/label/strg']):
            throats[id].label = strg.decode('utf-8')

    return network