Commit 4f3b6af2 authored by sfritschi's avatar sfritschi
Browse files

Further improved cell-list algo

parent da7acb9b
gen/plots/dendro.png

135 KB | W: | H:

gen/plots/dendro.png

82.2 KB | W: | H:

gen/plots/dendro.png
gen/plots/dendro.png
gen/plots/dendro.png
gen/plots/dendro.png
  • 2-up
  • Swipe
  • Onion skin
import multiprocessing as mp
def work(x, y):
res = []
for i in range(x):
res.append(y)
return res
def main():
with mp.Pool(processes=2) as p:
res = p.starmap(func=work, iterable=[(1, 2), (3, 4)])
print(res)
if __name__ == '__main__':
main()
......@@ -101,16 +101,41 @@ class Network:
return throat
class CellList:
def __init__(self, domainSize: List[float], cellSize: float, Lmax: float):
from math import floor
"""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, domainSize: List[float], cellSize: float, Lmax: float):
from math import ceil
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)]
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)
# Sort pores according to cell membership
self.poresSorted = [set() for _ in range(self.totalCells)]
for pore in pores:
cellIdx = self.pore_to_index(pore)
self.poresSorted[cellIdx].add(pore)
""" DEBUG
cellCounts = [0] * self.totalCells
for pore in pores:
cellIdx = self.pore_to_index(pore)
cellCounts[cellIdx] += 1
n = len(pores)
max_ = max(cellCounts)
min_ = min(cellCounts)
avg = sum(cellCounts) / len(cellCounts)
print("Cell statistics:")
print("Number of cells (%d, %d, %d)" % (self.nCells[0], self.nCells[1], self.nCells[2]))
print("Max. pores per cell: %d percent: %f" % (max_, max_ / n * 100.))
print("Min. pores per cell: %d percent: %f" % (min_, min_ / n * 100.))
print("Avg. pores per cell: %f percent: %f" % (avg, avg / n * 100.))
raise AssertionError("Stop")
"""
# Helper for converting 3D index triple into 1D (flattened) index
def flatten(self, i: int, j: int, k: int) -> int:
......@@ -129,14 +154,15 @@ class CellList:
for i in range(self.dim)])
return cellIdx
def valid_index(self, idx: int) -> bool:
def is_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 = []
# Compute dictionary of all nbor candidates of given pore
def find_candidates(self, pore: Pore) -> Dict[Pore, float]:
candidates = {}
cellIdxTrip = self.pore_to_triplet(pore)
# Check all possible neighbor cells (27 in 3D)
# TODO: Speed up finding of suitable candidates (e.g. multiprocessing)
for i in range(-1, 2):
for j in range(-1, 2):
for k in range(-1, 2):
......@@ -144,44 +170,20 @@ class CellList:
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 (Statistics)
"""
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 (Verification)
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
#assert(self.is_valid_index(nborIdx))
# Check all pores within current neighbor-cell
for nbor in self.poresSorted[nborIdx]:
if (nbor == pore): continue
l = distance(nbor.pos, pore.pos)
if (l < self.Lmax):
candidates[nbor] = l
return candidates
# Remove pore from sorted pores
def expel(self, pore: Pore):
poreIdx = self.pore_to_index(pore)
self.poresSorted[poreIdx].remove(pore)
def distance(p1: List[float], p2: List[float]):
"""Distance between points p1 and p2."""
......@@ -563,48 +565,43 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
"""
# 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)
# Initialize cell-list; Place each pore in respective cell-set
cellList = CellList(pores, trueDomainSize, basenet.Lmax, basenet.Lmax)
# throats, connect pores
if (not mute): print("\b"*21 + "connecting")
throats = []
k = 0 # count of unrealised throats
# Spawn parallel
for i, pore in enumerate(pores[:n]):
# (re)run triangulation if maximal number of neighbors exceeds
# threshold given by 2x init. maximum (large nnbrs -> expensive search)
# check throats and neighborhood
if (not mute):
print("\b"*24 + f"pore {i:7d}" + " "*12, end="", flush=True)
print(f"progress {((i+1) / n)*100.:.1f}%", end="\r", flush=True)
if (len(pore.throats) == 0): continue # no further connection needed
neighborhood = cellList.nbor_indices(pore)
conpores = set() # pores that pore is connected to
# find candidates in neighborhood and respective throat lengths
candidates = {}
for idx in neighborhood:
for nbor in poresSorted[idx]:
if (nbor == pore): continue
l = distance(pore.pos, nbor.pos)
if (l < basenet.Lmax):
candidates[nbor] = l
# TODO: Search neighboring cells in parallel
candidates = cellList.find_candidates(pore)
# establish throat connections between pore and best candidate
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()
lt = distance(*throat_ends(throat,
[ub-lb for lb,ub in zip(basenet.lb,basenet.ub)])) # target len.
# find best candidate
# avoid double connections and throats longer than Lmax
# TODO: Search neighboring cells in parallel
# TODO: Replace sort with linear search operation
ranking = [(l,p) for p, l in candidates.items() \
if (not p in conpores)]
if (len(ranking) == 0):
if (len(candidates) == 0):
k = k+1; continue # no candidates left, give up
lt, nbor = min(ranking, key=lambda lp: abs(lp[0] - lt))
nbor, lt = min(candidates.items(), key=lambda pl: abs(pl[1] - lt))
# Remove nbor from candidates
del candidates[nbor]
nborc = nbor # memorize potential periodic copy
# find original of buffer layer pore
if (nbor.label != LABELS[0]):
......@@ -614,33 +611,23 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
lbl = LABELS[0]
# connect pore to nbor
throats.append([pore, nbor, lbl, throat.r])
conpores.add(nborc) # update list of connected pores
# remove most similar throat of nbor
ranking = [(distance(*throat_ends(nthroat,
[ub-lb for lb,ub in zip(basenet.lb,basenet.ub)])),
nthroat) for nthroat in nbor.throats]
_, nbor_throat = min(ranking, key=lambda lth: abs(lth[0] - lt))
nbor.throats.remove(nbor_throat)
_, nthroat = min(ranking, key=lambda lth: abs(lth[0] - lt))
nbor.throats.remove(nthroat)
# remove fully connected nbor from search neighborhood
if (len(nbor.throats) == 0):
# remove fully-connected nbor pore from candidates
del candidates[nborc]
nborIdx = cellList.pore_to_index(nbor)
poresSorted[nborIdx].remove(nbor)
cellList.expel(nbor)
for cpore in copies[nbor]:
# Remove potential periodic copies of nbor as well
if cpore in candidates:
del candidates[cpore]
cnborIdx = cellList.pore_to_index(cpore)
poresSorted[cnborIdx].remove(cpore)
cellList.expel(cpore)
# remove fully connected pore from search neighborhood
poreIdx = cellList.pore_to_index(pore)
poresSorted[poreIdx].remove(pore)
cellList.expel(pore)
for cpore in copies[pore]:
cporeIdx = cellList.pore_to_index(cpore)
poresSorted[cporeIdx].remove(cpore)
cellList.expel(cpore)
print("\b"*24 + "{0:d} throats in total, {1:d} unrealised".\
format(len(throats), k))
......@@ -649,8 +636,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
network = Network(lb=[0.0 for k in range(d)],
ub=L, Lmax=basenet.Lmax, label='from_' + basenet.label)
for pore in pores[:n]: network.add_pore(pore)
for throat in throats: network.connect_pores(pore1=throat[0],
pore2=throat[1], label=throat[2], r=throat[3])
for throat in throats: network.connect_pores(*throat)
return network
......
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