Commit f7470925 authored by sfritschi's avatar sfritschi
Browse files

Matchmaker ignores pores already connected to current pore

parent 032a0f63
......@@ -27,7 +27,7 @@ def main():
start = perf_counter()
dendro = netflow.generate_dendrogram(basenet=basenet, targetsize=target, \
cutoff=cutoff, sd=42, nthreads=nthreads, \
maxIters=6, throatTol=0.01, mute=False)
mute=False)
end = perf_counter()
elapsed = end - start
print("Elapsed time: %e s" % elapsed)
......
......@@ -12,7 +12,7 @@ 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):
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
......@@ -23,6 +23,7 @@ class Pore:
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
......@@ -140,6 +141,8 @@ class CellList:
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)]
throatTotal = 0
for i, pore in enumerate(pores):
......@@ -200,18 +203,25 @@ class CellList:
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]:
if (nbor.index == pore.index): continue
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[nbor.index] = l
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[nbor.index] = l
candidates[nborIdx] = l
return candidates
......@@ -246,6 +256,10 @@ class CellList:
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]
......@@ -270,7 +284,11 @@ class CellList:
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
......@@ -635,7 +653,7 @@ def cell_list_scaling(basenet: Network, targetsize: List[int], \
def generate_dendrogram(basenet: Network, targetsize: List[int], \
cutoff: float = float('inf'), sd: int = None, nthreads: int = mp.cpu_count(), \
maxIters: int = 6, throatTol: float = 0.01, mute: bool = False) -> Network:
mute: bool = False) -> Network:
"""Generate a new spatially periodic network.
Based on an existing network basenet, generate a new network of size
......@@ -730,6 +748,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
# 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]
......@@ -739,7 +758,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
# 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
......@@ -748,21 +767,18 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
# Divide by two (since we are counting every throat twice)
totalThroats = len(cellList.poreMatchTable) // 2
# Evenly distribute pores (in interior) among threads
remainder = n % nthreads
loads = [n // nthreads] * nthreads
for i in range(remainder):
loads[i] += 1
# Compute inclusive-scan
displs = [0]
displs += accumulate(loads)
# 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.
......@@ -771,33 +787,48 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
# 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]]
# Set of realized throats Tuple[int, int]
realizedThroats = set()
if (not mute): print("Total throats: %d" % totalThroats)
iterCount = 0
throatsLeft = totalThroats
throatThreshold = int(throatTol * 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 (iterCount < maxIters and throatsLeft > throatThreshold):
while (throatsLeft != 0):
if (not mute):
print("Throats left: %d (%.1f%%)" % \
(throatsLeft, throatsLeft / totalThroats * 100.))
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)
print("Throats left: %d Throats Unrealized: %d (%.2f%%)" % \
(throatsLeft, throatsUnrealized, throatsUnrealized / 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)]
for worker in workers:
worker.join()
# 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):
......@@ -806,8 +837,9 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
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)
......@@ -827,23 +859,18 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
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]
# Make sure tuples of index pairs is unique
if (nbor.index < pore.index):
throatConn = (nbor.index, pore.index)
else:
throatConn = (pore.index, nbor.index)
# Check that nbor is not fully connected and that
# this pore is not already connected to it
if (len(nborTtbr) == 0 or throatConn in realizedThroats):
if (len(nborTtbr) == 0 or nbor.index in pConn):
continue
# Update set of realized throats
realizedThroats.add(throatConn)
# Update sets of connected pores
cellList.update_connection(pore.index, nbor.index)
r = Network.throatR[throatIdx]
targetLen = Network.throatL[throatIdx]
......@@ -880,46 +907,48 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
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
# Recompute loads and displacements for all threads
nremain = len(poresRemain)
remainder = nremain % nthreads
loads = [nremain // nthreads] * nthreads
for i in range(remainder):
loads[i] += 1
# Compute inclusive-scan
displs = [0]
displs += accumulate(loads)
throatsLeft = totalThroats - len(throats)
iterCount += 1
throatsLeft = throatCount // 2 # Counted all throats twice
throatsUnrealized = totalThroats - len(throats)
# Unrealized throats
nunrealized = throatsLeft
nUnrealized = throatsUnrealized
# Free memory (not needed anymore)
del cellList
del ttbrs
del realizedThroats
# DEBUG
if (iterCount == maxIters):
print("Maximum iteration count reached!")
percent = nunrealized / totalThroats * 100.
percent = nUnrealized / totalThroats * 100.
print("\b"*24 + "{0:d} throats in total, {1:d} unrealised ({2:.3f}%)".\
format(len(throats), nunrealized, percent))
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.))
......@@ -933,8 +962,8 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
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:
......@@ -956,7 +985,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
checkUnrealized //= 2 # counted every throat twice
# Make sure missing throats are accounted for
assert(checkUnrealized == nunrealized)
assert(checkUnrealized == nUnrealized)
"""
return network
......@@ -980,23 +1009,29 @@ def __add_buffer_layers(pores: List[Pore], targetsize: List[float],
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 = 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[int(lbl[0])]].add(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 = 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[int(lbl[0])]].add(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]
......
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