netgen.py 52.1 KB
Newer Older
1
2
3
4
5
# generation and planar cutting of periodic pore networks
#
# Daniel W. Meyer
# Institute of Fluid Dynamics, ETH Zurich
# January 2019
sfritschi's avatar
sfritschi committed
6
import multiprocessing as mp
7
8
9
10
11
12
13
14
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], \
15
        throats: Set = None, id: int = -1, index: int = -1):
16
17
18
19
20
21
22
23
24
        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'
sfritschi's avatar
sfritschi committed
25
26
        # Identifier for pores in list used in generate_dendrogram
        self.index = index
27
28
29
30
31
32
33
34
35
36
37
38
39
40
    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, \
41
        label: str = LABELS[0], id: int = -1, index: int = -1):
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
        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

105
class CellList:
sfritschi's avatar
sfritschi committed
106
107
108
109
    """CellList class dividing physical domain into equally-sized cells
    
    find_candidates() assumes pore is located in interior of domain (non-periodic)
    """
110
111
    def __init__(self, pores: List[Pore], copies: Dict[Pore, Set[Pore]],  \
                    domainSize: List[float], cellSize: float, \
112
                        basenet: Network, nInterior: int):
sfritschi's avatar
sfritschi committed
113
        from math import ceil
114
        self.dim = len(domainSize)
115
        self.Lmax = basenet.Lmax
116
        # Number of cells in each dimension
sfritschi's avatar
sfritschi committed
117
        self.nCells = [max(1, ceil(domainSize[i] / cellSize)) for i in range(self.dim)]
118
119
120
121
        # 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)
sfritschi's avatar
sfritschi committed
122
        # Pore positions 
123
124
125
126
127
128
129
130
131
        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
132
        # Used to keep track of already realized throats
133
        self.poreThroatIndices = [set() for _ in range(nInterior)]
134
            
135
        L = [ub-lb for lb,ub in zip(basenet.lb,basenet.ub)]
sfritschi's avatar
sfritschi committed
136
137
        # Sort pores according to cell membership
        self.poresSorted = [set() for _ in range(self.totalCells)]
sfritschi's avatar
sfritschi committed
138
        
139
140
141
142
143
144
145
146
147
        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)
148
                
149
150
                offset += nthroats
                self.throatOffsets.append(offset)
sfritschi's avatar
sfritschi committed
151
                    
152
153
154
155
                for throat in pore.throats:
                    self.throatL.append(distance(*throat_ends(throat, L)))
                    self.throatR.append(throat.r)
                    self.poreThroatIndices[poreIdx].add(throatIdx)
sfritschi's avatar
sfritschi committed
156
                    
157
                    throatIdx += 1
158
        
159
        throatTotal = offset
sfritschi's avatar
sfritschi committed
160
        self.poreMatchTable = mp.RawArray('i', throatTotal)
161
162
163
164
165
        
    # 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)
    
sfritschi's avatar
sfritschi committed
166
    def pore_to_triplet(self, poreIdx: int) -> List[int]:
167
168
        from math import floor
        # Shift pore position to be positive first (buffer layer)
sfritschi's avatar
sfritschi committed
169
170
        porePos = self.fetch_pos(poreIdx)
        return [floor( (porePos[i] + self.Lmax) * self.invCellSizes[i]) \
171
172
173
174
175
176
177
178
179
                        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
    
sfritschi's avatar
sfritschi committed
180
    def is_valid_index(self, idx: int) -> bool:
181
        return 0 <= idx < self.totalCells
182
183
184
185
    
    def match_maker(self, pores: List[int], tid: int):
        n = len(pores)
        for poreIdx in pores:
186
187
            if (tid == 0):
                print(f"progress {((poreIdx+1) / n)*100.:.1f}%", end="\r", flush=True)
sfritschi's avatar
sfritschi committed
188
            candidates = self.find_candidates(poreIdx)
189
190
191
192
            
            throats = self.fetch_throat_indices(poreIdx)
            for idx, throatIdx in enumerate(throats):
                lt = self.throatL[throatIdx]
193
                # Find best match among candidates
sfritschi's avatar
sfritschi committed
194
                matchIdx, _ = min(candidates.items(), key=lambda il: abs(il[1] - lt))
195
196
                # Remove match from candidates
                del candidates[matchIdx]
sfritschi's avatar
sfritschi committed
197
198
                # Put best match index in poreThroatTable
                self.set_throat(poreIdx, matchIdx, idx)
199
    
sfritschi's avatar
sfritschi committed
200
    # Compute dictionary of all nbor candidates of given pore
201
202
    def find_candidates(self, poreIdx: int) -> Dict[int, float]:    
        candidates  = {}
203
        
sfritschi's avatar
sfritschi committed
204
        porePos = self.fetch_pos(poreIdx)
sfritschi's avatar
sfritschi committed
205
206
        neighborhood = self.nbor_indices(poreIdx)
        nearestCellIdx = neighborhood[0]
207
        
sfritschi's avatar
sfritschi committed
208
        # Search current cell
sfritschi's avatar
sfritschi committed
209
210
211
212
        for nborIdx in self.poresSorted[nearestCellIdx]:
            if (nborIdx == poreIdx): continue
            l = distance(self.fetch_pos(nborIdx), porePos)
            if (l < self.Lmax):
213
                candidates[nborIdx] = l
sfritschi's avatar
sfritschi committed
214
215
216
217
218
219
        # 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
sfritschi's avatar
sfritschi committed
220
        
221
        return candidates
sfritschi's avatar
sfritschi committed
222
    
223
224
225
    def set_throat(self, poreIdx: int, matchIdx: int, idx: int):
        index = self.throatOffsets[poreIdx] + idx
        # Set matchIdx
sfritschi's avatar
sfritschi committed
226
        self.poreMatchTable[index] = matchIdx
227
        
sfritschi's avatar
sfritschi committed
228
    def fetch_matches(self, poreIdx: int) -> List[int]:
229
230
        lb = self.throatOffsets[poreIdx]
        ub = self.throatOffsets[poreIdx + 1]
sfritschi's avatar
sfritschi committed
231
        return self.poreMatchTable[lb:ub]
232
233
234
235
236
237
    
    def fetch_throat_indices(self, poreIdx: int) -> List[int]:
        lb = self.throatOffsets[poreIdx]
        ub = self.throatOffsets[poreIdx + 1]
        return range(lb, ub)
        
sfritschi's avatar
sfritschi committed
238
239
240
241
242
    def fetch_pos(self, poreIdx: int) -> List[float]:
        lb = self.dim * poreIdx
        ub = self.dim + lb
        return self.poresPos[lb:ub]
    
243
    # Consider all 27 neighboring cells
244
245
    def nbor_indices(self, poreIdx: int) -> List[int]:
        cellIdx = self.pore_to_triplet(poreIdx)
sfritschi's avatar
sfritschi committed
246
247
248
        base   = self.flatten(*cellIdx)
        widthx = self.nCells[0]
        widthy = widthx * self.nCells[1]
249
250
251
252
253
254
255
256
257
258
259
260
261
        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]
sfritschi's avatar
sfritschi committed
262
263
264
265
        
    def nbor_cell_counts(self, neighborhood: List[int]) -> List[int]:
        return [len(self.poresSorted[i]) for i in neighborhood]
        
sfritschi's avatar
sfritschi committed
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
    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)

sfritschi's avatar
sfritschi committed
282
283
    # Remove pore from sorted pores
    def expel(self, pore: Pore):
sfritschi's avatar
sfritschi committed
284
285
        poreCellIdx = self.pore_to_index(pore)
        self.poresSorted[poreCellIdx].remove(pore.index)
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555

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
sfritschi's avatar
sfritschi committed
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647

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)
    
648
    
649
def generate_dendrogram(basenet: Network, targetsize: List[int], \
650
651
    cutoff: float = float('inf'), sd: int = None, nthreads: int = mp.cpu_count(), \
    mute: bool = False) -> Network:
652
653
654
655
656
657
658
659
660
661
662
    """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.
    """
sfritschi's avatar
sfritschi committed
663
664
    import multiprocessing as mp
    from random import seed, random, randint
665
    from itertools import product, accumulate
666
667
668
669
670
671
672
673
674
675
676
677
    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
sfritschi's avatar
sfritschi committed
678
    if (not mute): print("distributing pores...")
679
680
681
682
683
684
685
    # 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
686
        print("uniform pore distribution")
687
        # loop over network sections & add uniformly distributed pores drawn from basenet
sfritschi's avatar
sfritschi committed
688
        runningIndex = 0
689
690
691
692
        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], 
sfritschi's avatar
sfritschi committed
693
694
                    throats=pore.throats.copy(), index=runningIndex))
                runningIndex += 1
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
    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)]
sfritschi's avatar
sfritschi committed
731
732
                pores.insert(randint(0,len(pores)), Pore(pos=pos,
                    r=pore.r, label=LABELS[0], throats=pore.throats.copy()))
sfritschi's avatar
sfritschi committed
733
                
734
        print("left {0:d} of {1:d} clusters (incl. {2:d} pores) untouched".\
735
            format(sum([int(not j) for j in touched]), len(touched), len(basepores)))
sfritschi's avatar
sfritschi committed
736
737
738
        # flip pores outside back into domain and set respective index of all pores
        for i, pore in enumerate(pores):
            pore.index = i
739
740
741
742
743
744
745
746
747
            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)
748
    
749
    isPeriodicCheck = 1 in targetsize
750
751
    # Domain size including periodic buffer layers on all sides
    trueDomainSize = [L[i] + 2 * basenet.Lmax for i in range(d)]
sfritschi's avatar
sfritschi committed
752
    # Initialize cell-list; Place each pore in respective cell-set
753
    cellList = CellList(pores, copies, trueDomainSize, basenet.Lmax, \
754
                            basenet, n)
sfritschi's avatar
sfritschi committed
755
    
756
    # Evenly distribute pores (in interior) among threads
757
758
759
760
    remainder = n % nthreads
    loads = [n // nthreads] * nthreads
    for i in range(remainder):
        loads[i] += 1
sfritschi's avatar
sfritschi committed
761
    
762
763
764
    # Compute inclusive-scan
    displs = [0]
    displs += accumulate(loads)
765
766
    
    workers = []
767
    if (not mute): print("computing best matches")
768
    # Compute best matches for each pore in parallel using available threads
769
770
771
    for tid in range(nthreads):
        poreList = range(displs[tid], displs[tid+1])
        worker = mp.Process(target=cellList.match_maker, args=(poreList, tid))
sfritschi's avatar
sfritschi committed
772
        worker.start()
773
        workers.append(worker)
774
775
776
    
    for worker in workers:
        worker.join()
777
        
778
    # throats, connect pores
sfritschi's avatar
sfritschi committed
779
    if (not mute): print("\nconnecting")
780
781
782
    throats = []
    k = 0 # count of unrealised throats
    
783
784
    avg_throat_diff = 0.
    count = 0
785
786
    total = 0
    n_already_taken = 0
787
    for poreIdx, pore in enumerate(pores[:n]):
788
789
        # Delete throats (not needed anymore)
        pore.throats.clear()
790
        
791
        if (not mute):
792
            print(f"progress {((poreIdx+1) / n)*100.:.1f}%", end="\r", flush=True)
793
794
795
        # Set of throats to be realized
        ttbr = cellList.poreThroatIndices[poreIdx]
        if (len(ttbr) == 0): continue  # Nothing left to do
sfritschi's avatar
sfritschi committed
796
        
sfritschi's avatar
sfritschi committed
797
        matches = cellList.fetch_matches(poreIdx)
798
        throatIndices = cellList.fetch_throat_indices(poreIdx)
799
        
800
        conporesInd = set()  # Set of connected pores (needed if target size contains 1)
801
802
        for i, match in enumerate(matches):
            throatIdx = throatIndices[i]
803
            
804
805
            if (throatIdx not in ttbr):
                continue  # Throat already realized
sfritschi's avatar
sfritschi committed
806
                            
807
            nbor = pores[match]
sfritschi's avatar
sfritschi committed
808
            
809
810
811
812
813
814
815
            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]
sfritschi's avatar
sfritschi committed
816
            
817
            nborTtbr = cellList.poreThroatIndices[nbor.index]
818
            
819
            total += 1
820
            if (len(nborTtbr) == 0 or nbor.index in conporesInd):
821
822
823
                n_already_taken += 1
                continue
            
824
825
            conporesInd.add(nbor.index)
            
826
827
828
829
830
            r = cellList.throatR[throatIdx]
            targetLen = cellList.throatL[throatIdx]
            # Compute actual throat length between pores
            lt = distance(pore.pos, nborc.pos)
            
831
832
833
834
835
836
837
838
            # 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)
839
840
            avg_throat_diff += diff
            count += 1
sfritschi's avatar
sfritschi committed
841
            
842
843
844
845
846
            # Remove throat from neighbor
            nborTtbr.remove(nborThroatIdx)
            # Successfully connected this throat
            ttbr.remove(throatIdx)
            
847
            # remove fully-connected nbor from search neighborhood
848
849
850
851
852
853
            if (len(nborTtbr) == 0):
                cellList.expel(nbor)
                
                for cpore in copies[nbor]:
                    cellList.expel(cpore)
        
854
        # remove fully-connected pore from search neighborhood
855
856
857
858
859
        if (len(ttbr) == 0):
            cellList.expel(pore)
        
            for cpore in copies[pore]:
                cellList.expel(cpore)
sfritschi's avatar
sfritschi committed
860
861
    
    # Free space
sfritschi's avatar
sfritschi committed
862
    del cellList.poreMatchTable
sfritschi's avatar
sfritschi committed
863
    
864
    # Finally, find next best candidates for left-out pores
865
    # in REVERSED order
sfritschi's avatar
sfritschi committed
866
    if (not mute): print("\nfinding next best candidates")
867
    
868
869
870
    poreIdx = n
    for it, pore in enumerate(reversed(pores[:n])):
        poreIdx -= 1
871
872
        
        if (not mute):
873
            print(f"progress {((it+1) / n)*100.:.1f}%", end="\r", flush=True)
874
875
876
877
        
        ttbr = cellList.poreThroatIndices[poreIdx]
        if (len(ttbr) == 0):
            continue  # Nothing left to do
sfritschi's avatar
sfritschi committed
878
879
        
        # Expensive: Find alternative (multiprocessing)
880
881
        candidates = cellList.find_candidates(poreIdx)
            
sfritschi's avatar
sfritschi committed
882
        for throatIdx in ttbr.copy():
883
            
sfritschi's avatar
sfritschi committed
884
            r = cellList.throatR[throatIdx]
885
886
887
888
889
890
891
892
            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))
893
                        
894
            nbor = pores[nborIdx]
895
896
897
898
899
900
901
            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]
sfritschi's avatar
sfritschi committed
902
            
903
904
905
906
907
908
909
910
911
912
913
914
            # 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]

915
            # connect pore to nbor
916
917
918
            throats.append([pore, nbor, lbl, r])
            
            nborTtbr = cellList.poreThroatIndices[nbor.index]
919
            # remove most similar throat of nbor
920
921
922
            ranking = [(nborThroatIdx, abs(cellList.throatL[nborThroatIdx] - lt)) \
                         for nborThroatIdx in nborTtbr]
            nborThroatIdx, _ = min(ranking, key=lambda il: il[1])
923
            
924
925
926
927
            diff = abs(targetLen - lt)
            avg_throat_diff += diff
            count += 1
            
sfritschi's avatar
sfritschi committed
928
929
            # Remove throat from pore
            ttbr.remove(throatIdx)
930
931
            # Remove throat from neighbor
            nborTtbr.remove(nborThroatIdx)
932
            
933
            # remove fully connected nbor from search neighborhood
934
            if (len(nborTtbr) == 0):
sfritschi's avatar
sfritschi committed
935
                cellList.expel(nbor)
sfritschi's avatar
sfritschi committed
936
                
937
                for cpore in copies[nbor]:
sfritschi's avatar
sfritschi committed
938
                    cellList.expel(cpore)
939
        
940
        # remove fully connected pore from search neighborhood
941
942
        if (len(ttbr) == 0):
            cellList.expel(pore)
943
        
944
945
            for cpore in copies[pore]:
                cellList.expel(cpore)
946
    
947
948
949
    # Free memory associated with cellList (not needed anymore)
    del cellList
    
sfritschi's avatar
sfritschi committed
950
    # DEBUG
sfritschi's avatar
sfritschi committed
951
952
953
    percent = k / (k + len(throats)) * 100.
    print("\b"*24 + "{0:d} throats in total, {1:d} unrealised ({2:.1f}%)".\
        format(len(throats), k, percent))
954
955
956
    avg_throat_diff /= count
    print("Avg. throat length difference: %e" % avg_throat_diff)
    print("Relative to Lmax: %.1f%%" % (avg_throat_diff / basenet.Lmax * 100.))
957
    print("Percentage where match was fully-connected: %.1f%%" % (n_already_taken / total * 100.))
958

959
960
961
    # assemble and return network
    network = Network(lb=[0.0 for k in range(d)],
        ub=L, Lmax=basenet.Lmax, label='from_' + basenet.label)
962
963
    for pore in pores[:n]:
        network.add_pore(pore)
964
965
    for throat in throats: network.connect_pores(pore1=throat[0],
        pore2=throat[1], label=throat[2], r=throat[3])
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
    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}}
sfritschi's avatar
sfritschi committed
982
    runningIndex = n
983
984
985
986
987
    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:
sfritschi's avatar
sfritschi committed
988
            pcopy = Pore(pore.pos, pore.r, throats=pore.throats, index=runningIndex, id=pore.id)
989
990
991
992
993
            pcopy.pos[k] = pcopy.pos[k]-targetsize[k]
            lbl = pore.label.split(); lbl[k + 1] = '-1'
            pcopy.label = ' '.join(lbl) # periodicity label
            pores.append(pcopy)
            copies[pores[int(lbl[0])]].add(pcopy)
sfritschi's avatar
sfritschi committed
994
            runningIndex += 1
995
996
997
998
        # right bound
        pores_layer = [pore for pore in pores \
            if 0 <= pore.pos[k] < Lbuffer]
        for pore in pores_layer:
sfritschi's avatar
sfritschi committed
999
            pcopy = Pore(pore.pos, pore.r, throats=pore.throats, index=runningIndex, id=pore.id)
1000
            pcopy.pos[k] = pcopy.pos[k]+targetsize[k]