Commit d933ecdc authored by sfritschi's avatar sfritschi
Browse files

Replaced triangulation with cell lists

parent 9b2b4264
......@@ -67,12 +67,12 @@ def main():
assert(len(pore.throats) == 0)
"""
target = [1, 1, 1]
target = [3, 3, 3]
cutoff = 0.5 * max([basenet.ub[i] - basenet.lb[i] \
for i in range(len(basenet.ub))])
dendro = netflow.generate_dendrogram(basenet=basenet, targetsize=target, \
cutoff=cutoff, sd=42, mute=True)
cutoff=cutoff, sd=42, mute=False)
_ = netgen.plot_network(dendro)
plt.savefig('plots/dendro.png')
......
gen/plots/dendro.png

142 KB | W: | H:

gen/plots/dendro.png

135 KB | W: | H:

gen/plots/dendro.png
gen/plots/dendro.png
gen/plots/dendro.png
gen/plots/dendro.png
  • 2-up
  • Swipe
  • Onion skin
......@@ -8,8 +8,6 @@ from typing import List, Dict, Set, Tuple # for type hints in argument lists
LABELS = ('', 'in', 'out', 'cut') # don't change the order
import multiprocessing as mp
class Pore:
"""Class of a pore connected to other pores via throats."""
id = 0
......@@ -24,7 +22,6 @@ class Pore:
self.throats = throats # throats set
self.r = r # radius
self.label = label # label string like '', 'in', or 'out'
#self.is_connected = mp.Value('i', 0)
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
......@@ -374,7 +371,89 @@ def generate_imperial(basenet: Network, targetsize: List[float], \
pore2=pores[throat[1]], label=throat[2], r=throat[3])
return network
class CellList:
def __init__(self, domainSize: List[float], cellSize: float, Lmax: float):
from math import floor
self.dim = len(domainSize)
self.Lmax = Lmax
# Number of cells in each dimension
self.nCells = [max(1, floor(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)
# 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) -> List[int]:
from math import floor
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
def valid_index(self, idx: int) -> bool:
return 0 <= idx < self.totalCells
# Compute list of indices of (valid) nbor pores
def nbor_indices(self, pore: Pore) -> List[int]:
nborIndices = []
cellIdxTrip = self.pore_to_triplet(pore)
# Check all possible neighbor cells (27 in 3D)
for i in range(-1, 2):
for j in range(-1, 2):
for k in range(-1, 2):
current = [cellIdxTrip[0] + i, cellIdxTrip[1] + j, \
cellIdxTrip[2] + k]
nborIdx = self.flatten(*current)
# Check if valid cell index
if self.valid_index(nborIdx):
nborIndices.append(nborIdx)
return nborIndices
# Compute cell index of each pore and add it to respective set
def sort_pores(self, pores: Set[Pore]) -> List[Set]:
# DEBUG
"""
cellCounts = [0] * self.totalCells
for pore in pores:
cellIdx = self.pore_to_index(pore)
cellCounts[cellIdx] += 1
n = len(pores)
max_ = max(cellCounts)
min_ = min(cellCounts)
avg = sum(cellCounts) / len(cellCounts)
print("Cell statistics:")
print("Number of cells (%d, %d, %d)" % (self.nCells[0], self.nCells[1], self.nCells[2]))
print("Max. pores per cell: %d percent: %f" % (max_, max_ / n * 100.))
print("Min. pores per cell: %d percent: %f" % (min_, min_ / n * 100.))
print("Avg. pores per cell: %f percent: %f" % (avg, avg / n * 100.))
"""
#raise AssertionError("Stop")
poresSorted = [set() for _ in range(self.totalCells)]
for pore in pores:
cellIdx = self.pore_to_index(pore)
poresSorted[cellIdx].add(pore)
"""
# DEBUG
for i in range(self.totalCells):
for pore in poresSorted[i]:
cellIdx = self.pore_to_index(pore)
if cellIdx != i:
print((cellIdx, i))
raise AssertionError("Wrong index")
"""
return poresSorted
def generate_dendrogram(basenet: Network, targetsize: List[int], \
cutoff: float = float('inf'), sd: int = None, mute: bool = False) \
-> Network:
......@@ -509,80 +588,34 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
nnbrs = max(nnbrs, len(neighborhood[nbor])) # update max. number of neighbors
del(neighborhood[outsider])
return nnbrs
# Worker/consumer function
def process_throat(input_queue, results_queue):
# Infinite work-loop
while True:
# Receive throat (blocking)
input_ = input_queue.get()
throat = input_[0]
front = input_[1]
neighborhood = input_[2]
candidates = input_[3]
lt = distance(*throat_ends(throat,
[ub-lb for lb,ub in zip(basenet.lb,basenet.ub)])) # target len.
# gather candidates
while (len(front) > 0):
if (min(front.values()) >= lt): break
# grow search region with candidates
newfront = {}
for frontpore in front:
for nbor in neighborhood[frontpore]:
if (nbor == pore): continue # no self connections
l = distance(pore.pos, nbor.pos)
# distinguish between new front-/candidate-pores
if (not nbor in candidates): newfront[nbor] = l
candidates[nbor] = l
front = {p:l for p, l in newfront.items()} # copy
# find best candidate
# avoid double connections and throats longer than Lmax
ranking = [(l,p) for p, l in candidates.items() \
if (True) and (l < basenet.Lmax)]
if (len(ranking) == 0):
continue # no candidates left, give up
ranking.sort(key=lambda lp: abs(lp[0]-lt))
nbor = ranking[0][1] # pore
nborc = nbor # memorize potential periodic copy
# find original of buffer layer pore
if (nbor.label != LABELS[0]):
j, lbl = nbor.label.split(' ',1)
nbor = pores[int(j)]
else:
lbl = LABELS[0]
# remove most similar throat of nbor
lt = ranking[0][0] # actual throat length
ranking = [(distance(*throat_ends(nthroat,
[ub-lb for lb,ub in zip(basenet.lb,basenet.ub)])),
nthroat) for nthroat in nbor.throats]
ranking.sort(key=lambda lth: abs(lth[0]-lt))
# Put results in queue
results_queue.put([nbor, lbl, throat.r, ranking[0][1]])
"""
DEBUG: CELL LIST STATISTICS
print(len(pores))
# Optimal cutoff; no interaction beyond cutoff
cutoff = max([distance(*throat_ends(throat,
[ub-lb for lb,ub in zip(basenet.lb,basenet.ub)])) for pore in pores for throat in pore.throats])
print("Cutoff: %e, Lmax: %e" % (cutoff, basenet.Lmax))
from collections import Counter
cellmembers = Counter([(int(pore.pos[0] / cutoff), int(pore.pos[1] / cutoff), int(pore.pos[2] / cutoff)) for pore in pores])
print("Max: %d" % max(cellmembers.values()))
print("Min: %d" % min(cellmembers.values()))
print("Avg: %f" % (sum(cellmembers.values()) / len(cellmembers)))
"""
# Domain size including periodic buffer layers on all sides
trueDomainSize = [L[i] + 2 * basenet.Lmax for i in range(d)]
cellList = CellList(trueDomainSize, basenet.Lmax, basenet.Lmax)
poresSorted = cellList.sort_pores(pores)
# throats, connect pores
if (not mute): print("\b"*21 + "connecting")
throats = []
n_unrealized = 0 # count of unrealised throats
nnbrs0 = 0; nnbrs = 0; # (initial) maximum number of neighbors
# Spawn processes to process individual throats
import sys
# Required for pickling recursive Throat object
sys.setrecursionlimit(8000)
import multiprocessing as mp
input_queue = mp.Queue() # Queue for sending throats to workers
results_queue = mp.Queue() # Queue for sending results back
nworkers = 4
for _ in range(nworkers):
proc = mp.Process(target=process_throat, args=(input_queue, results_queue), daemon=True)
proc.start()
k = 0 # count of unrealised throats
#nnbrs0 = 0; nnbrs = 0; # (initial) maximum number of neighbors
for i, pore in enumerate(pores[:n]):
print("Processing pore: %d/%d" % (i, n+1), end="\r", flush=True)
# (re)run triangulation if maximal number of neighbors exceeds
# threshold given by 2x init. maximum (large nnbrs -> expensive search)
"""
if ((i == 0) or (nnbrs > 2*nnbrs0)):
if (not mute): print(" triangulation", end="", flush=True)
import scipy.spatial.qhull as qh
......@@ -593,10 +626,12 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
print("Error: triangulation failed,", err)
if (i == 0): raise
if (not mute): print("\b"*14 + " "*14 + "\b"*14, end="", flush=True)
"""
# check throats and neighborhood
if (not mute):
print("\b"*24 + f"pore {i:7d}" + " "*12, end="", flush=True)
if (len(pore.throats) == 0): continue # no further connection needed
"""
if neighborhood == {}: break # no pores left to connect to
if (i == 0): nnbrs0 = nnbrs
# lists with front and candidate pores
......@@ -604,40 +639,88 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
front = dict(zip(front,
[distance(pore.pos, nbor.pos) for nbor in front]))
candidates = {p:l for p, l in front.items()} # copy front dict
#conpores = set() # pores that pore is connected to
"""
neighborhood = cellList.nbor_indices(pore)
conpores = set() # pores that pore is connected to
# establish throat connections between pore and best candidate
nthroats = len(pore.throats)
while (len(pore.throats) > 0):
if (not mute):
print("\b"*12 + ", throat {0:3d}".format(len(pore.throats)),
end="", flush=True)
throat = pore.throats.pop()
# Put throat in queue to be processed by workers
input_queue.put([throat, front, neighborhood, candidates])
for _ in range(nthroats):
# Fetch results (blocking)
results = results_queue.get()
nbor = results[0]
lbl = results[1]
throat_r = results[2]
nbor_throat = results[3]
lt = distance(*throat_ends(throat,
[ub-lb for lb,ub in zip(basenet.lb,basenet.ub)])) # target len.
"""
# gather candidates
while (len(front) > 0):
if (min(front.values()) >= lt): break
# grow search region with candidates
newfront = {}
for frontpore in front:
for nbor in neighborhood[frontpore]:
if (nbor == pore): continue # no self connections
l = distance(pore.pos, nbor.pos)
# distinguish between new front-/candidate-pores
if (not nbor in candidates): newfront[nbor] = l
candidates[nbor] = l
front = {p:l for p, l in newfront.items()} # copy
"""
# find best candidate
# avoid double connections and throats longer than Lmax
candidates = {}
for idx in neighborhood:
for nbor in poresSorted[idx]:
if (nbor == pore): continue
l = distance(pore.pos, nbor.pos)
candidates[nbor] = l
# PROBLEM: pore in poresSorted should not have 0 throats -> inconsistent state
ranking = [(l,p) for p, l in candidates.items() \
if (not p in conpores) and (l < basenet.Lmax)]
if (len(ranking) == 0):
k = k+1; continue # no candidates left, give up
ranking.sort(key=lambda lp: abs(lp[0]-lt))
nbor = ranking[0][1] # pore
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]
# connect pore to nbor
throats.append([pore, nbor, lbl, throat_r])
#conpores.add(nborc) # update list of connected pores
nbor.throats.remove(nbor_throat)
throats.append([pore, nbor, lbl, throat.r])
conpores.add(nborc) # update list of connected pores
# remove most similar throat of nbor
lt = ranking[0][0] # actual throat length
ranking = [(distance(*throat_ends(nthroat,
[ub-lb for lb,ub in zip(basenet.lb,basenet.ub)])),
nthroat) for nthroat in nbor.throats]
ranking.sort(key=lambda lth: abs(lth[0]-lt))
nbor.throats.remove(ranking[0][1])
# remove fully connected nbor from search neighborhood
if (len(nbor.throats) == 0):
nnbrs = expel(pore, nbor, neighborhood, front, candidates, nnbrs)
nborIdx = cellList.pore_to_index(nbor)
poresSorted[nborIdx].remove(nbor)
#nnbrs = expel(pore, nbor, neighborhood, front, candidates, nnbrs)
for cpore in copies[nbor]:
nnbrs = expel(pore, cpore, neighborhood, front, candidates, nnbrs)
cnborIdx = cellList.pore_to_index(cpore)
poresSorted[cnborIdx].remove(cpore)
#nnbrs = expel(pore, cpore, neighborhood, front, candidates, nnbrs)
# remove fully connected pore from search neighborhood
nnbrs = expel(pore, pore, neighborhood, front, candidates, nnbrs)
poreIdx = cellList.pore_to_index(pore)
poresSorted[poreIdx].remove(pore)
#nnbrs = expel(pore, pore, neighborhood, front, candidates, nnbrs)
for cpore in copies[pore]:
nnbrs = expel(pore, cpore, neighborhood, front, candidates, nnbrs)
cporeIdx = cellList.pore_to_index(cpore)
poresSorted[cporeIdx].remove(cpore)
#nnbrs = expel(pore, cpore, neighborhood, front, candidates, nnbrs)
print("\b"*24 + "{0:d} throats in total, {1:d} unrealised".\
format(len(throats), n_unrealized))
format(len(throats), k))
# assemble and return network
network = Network(lb=[0.0 for k in range(d)],
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment