Commit 0e2123a4 authored by sfritschi's avatar sfritschi
Browse files

Fixed implementation using throats of pores directly by copy

parent 683437c3
......@@ -12,7 +12,7 @@ def main():
nthreads = int(sys.argv[1])
print("Using %d threads" % nthreads)
basenet = netflow.load_network_from('../netflow/network/network.h5')
basenet = netflow.load_network_from('../netflow/network/network.h5', isGenerator=True)
print("Network statistics:")
print("Throats: %d" % len(basenet.throats))
pore_throat_counts = [len(pore.throats) for pore in basenet.pores]
......
......@@ -22,7 +22,6 @@ class Pore:
self.throats = throats # throats set
self.r = r # radius
self.label = label # label string like '', 'in', or 'out'
# Identifier for pores in list used in generate_dendrogram
self.index = index
def __repr__(self): return str(self.__class__) + ': ' + \
str({k:self.__dict__[k] for k in self.__dict__ if k != 'throats'}) + \
......@@ -56,6 +55,10 @@ class Network:
Lmax must be smaller than the smallest side length of the network ub-lb.
"""
throatIdx = 0
# Globals
throatL = []
throatR = []
def __init__(self, lb: List[float] = [], ub: List[float] = [], \
pores: List[Pore] = None, throats: List[Throat] = None, \
Lmax: float = 0.0, label: str = None):
......@@ -73,6 +76,8 @@ class Network:
self.label = datetime.today().isoformat(sep=' ',timespec='seconds')
else:
self.label = label
# Lengths of network domain
self.L = [ub-lb for lb,ub in zip(self.lb,self.ub)]
def __repr__(self):
# count in/out resp. labelled pores
k = sum([int(pore.label != LABELS[0]) for pore in self.pores])
......@@ -101,14 +106,25 @@ class Network:
self.throats.add(throat)
pore1.throats.add(throat); pore2.throats.add(throat)
return throat
def connect_pores_generator(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)
lt = distance(*throat_ends(throat, self.L))
Network.throatL.append(lt)
Network.throatR.append(r)
self.throats.add(throat)
pore1.throats.add(Network.throatIdx)
pore2.throats.add(Network.throatIdx)
Network.throatIdx += 1
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], copies: Dict[Pore, Set[Pore]], \
domainSize: List[float], cellSize: float, \
def __init__(self, pores: List[Pore], domainSize: List[float], cellSize: float, \
basenet: Network, nInterior: int):
from math import ceil
self.dim = len(domainSize)
......@@ -119,44 +135,20 @@ class CellList:
self.invCellSizes = [self.nCells[i] / domainSize[i] for i in range(self.dim)]
# Total number of cells
self.totalCells = prod(self.nCells)
# Pore positions
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
# Used to keep track of already realized throats
self.poreThroatIndices = [set() for _ in range(nInterior)]
L = [ub-lb for lb,ub in zip(basenet.lb,basenet.ub)]
# Offsets in throat arrays
self.throatOffsets = [0]
# Sort pores according to cell membership
self.poresSorted = [set() for _ in range(self.totalCells)]
for poreIdx, pore in enumerate(pores):
throatTotal = 0
for i, pore in enumerate(pores):
cellIdx = self.pore_to_index(pore)
self.poresSorted[cellIdx].add(pore)
self.poresSorted[cellIdx].add(poreIdx)
self.poresPos += pore.pos
if poreIdx < nInterior:
# Pre-fix sum
nthroats = len(pore.throats)
if i < nInterior:
throatTotal += len(pore.throats)
self.throatOffsets.append(throatTotal)
offset += nthroats
self.throatOffsets.append(offset)
for throat in pore.throats:
self.throatL.append(distance(*throat_ends(throat, L)))
self.throatR.append(throat.r)
self.poreThroatIndices[poreIdx].add(throatIdx)
throatIdx += 1
throatTotal = offset
self.poreMatchTable = mp.RawArray('i', throatTotal)
# Helper for converting 3D index triple into 1D (flattened) index
......@@ -180,43 +172,42 @@ class CellList:
def is_valid_index(self, idx: int) -> bool:
return 0 <= idx < self.totalCells
def match_maker(self, pores: List[int], tid: int):
def match_maker(self, pores: List[Pore], tid: int):
n = len(pores)
for poreIdx in pores:
for i, pore in enumerate(pores):
if (tid == 0):
print(f"progress {((poreIdx+1) / n)*100.:.1f}%", end="\r", flush=True)
candidates = self.find_candidates(poreIdx)
print(f"progress {((i+1) / n)*100.:.1f}%", end="\r", flush=True)
candidates = self.find_candidates(pore)
throats = self.fetch_throat_indices(poreIdx)
for idx, throatIdx in enumerate(throats):
lt = self.throatL[throatIdx]
for idx, throatIdx in enumerate(pore.throats):
lt = Network.throatL[throatIdx]
# 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 poreThroatTable
self.set_throat(poreIdx, matchIdx, idx)
self.set_throat(pore.index, matchIdx, idx)
# Compute dictionary of all nbor candidates of given pore
def find_candidates(self, poreIdx: int) -> Dict[int, float]:
def find_candidates(self, pore: Pore) -> Dict[int, float]:
candidates = {}
porePos = self.fetch_pos(poreIdx)
neighborhood = self.nbor_indices(poreIdx)
porePos = pore.pos
neighborhood = self.nbor_indices(pore)
nearestCellIdx = neighborhood[0]
# Search current cell
for nborIdx in self.poresSorted[nearestCellIdx]:
if (nborIdx == poreIdx): continue
l = distance(self.fetch_pos(nborIdx), porePos)
for nbor in self.poresSorted[nearestCellIdx]:
if (nbor.index == pore.index): continue
l = distance(nbor.pos, porePos)
if (l < self.Lmax):
candidates[nborIdx] = l
candidates[nbor.index] = l
# Search neighboring cells
for nborCellIdx in neighborhood[1:]:
for nborIdx in self.poresSorted[nborCellIdx]:
l = distance(self.fetch_pos(nborIdx), porePos)
for nbor in self.poresSorted[nborCellIdx]:
l = distance(nbor.pos, porePos)
if (l < self.Lmax):
candidates[nborIdx] = l
candidates[nbor.index] = l
return candidates
......@@ -241,9 +232,8 @@ class CellList:
return self.poresPos[lb:ub]
# Consider all 27 neighboring cells
def nbor_indices(self, poreIdx: int) -> List[int]:
cellIdx = self.pore_to_triplet(poreIdx)
base = self.flatten(*cellIdx)
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,
......@@ -282,7 +272,7 @@ class CellList:
# Remove pore from sorted pores
def expel(self, pore: Pore):
poreCellIdx = self.pore_to_index(pore)
self.poresSorted[poreCellIdx].remove(pore.index)
self.poresSorted[poreCellIdx].remove(pore)
def distance(p1: List[float], p2: List[float]):
"""Distance between points p1 and p2."""
......@@ -683,9 +673,9 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
# 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
runningIndex = 0
for s in product(*[range(k) for k in targetsize]):
for k, pore in enumerate(basepores):
pos = [random()*l for l in L]
......@@ -745,13 +735,12 @@ 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)
isPeriodicCheck = 1 in targetsize
# 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, copies, trueDomainSize, basenet.Lmax, \
basenet, n)
cellList = CellList(pores, trueDomainSize, basenet.Lmax, basenet, n)
# Evenly distribute pores (in interior) among threads
remainder = n % nthreads
......@@ -767,14 +756,14 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
if (not mute): print("computing best matches")
# Compute best matches for each pore in parallel using available threads
for tid in range(nthreads):
poreList = range(displs[tid], displs[tid+1])
poreList = pores[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()
# throats, connect pores
if (not mute): print("\nconnecting")
throats = []
......@@ -784,26 +773,24 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
count = 0
total = 0
n_already_taken = 0
ttbrs = [pore.throats.copy() for pore in pores]
for poreIdx, pore in enumerate(pores[:n]):
# Delete throats (not needed anymore)
pore.throats.clear()
if (not mute):
print(f"progress {((poreIdx+1) / n)*100.:.1f}%", end="\r", flush=True)
# Set of throats to be realized
ttbr = cellList.poreThroatIndices[poreIdx]
if (len(ttbr) == 0): continue # Nothing left to do
poreTtbr = ttbrs[poreIdx]
if (len(poreTtbr) == 0): continue # Nothing left to do
matches = cellList.fetch_matches(poreIdx)
throatIndices = cellList.fetch_throat_indices(poreIdx)
conporesInd = set() # Set of connected pores (needed if target size contains 1)
for i, match in enumerate(matches):
throatIdx = throatIndices[i]
for throatIdx, match in zip(pore.throats, matches):
if (throatIdx not in ttbr):
if throatIdx not in poreTtbr:
continue # Throat already realized
nbor = pores[match]
nborc = nbor # memorize potential periodic copy
......@@ -814,7 +801,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
else:
lbl = LABELS[0]
nborTtbr = cellList.poreThroatIndices[nbor.index]
nborTtbr = ttbrs[nbor.index]
total += 1
if (len(nborTtbr) == 0 or nbor.index in conporesInd):
......@@ -823,15 +810,15 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
conporesInd.add(nbor.index)
r = cellList.throatR[throatIdx]
targetLen = cellList.throatL[throatIdx]
r = Network.throatR[throatIdx]
targetLen = Network.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(cellList.throatL[nborThroatIdx] - lt)) \
ranking = [(nborThroatIdx, abs(Network.throatL[nborThroatIdx] - lt)) \
for nborThroatIdx in nborTtbr]
nborThroatIdx, _ = min(ranking, key=lambda il: il[1])
......@@ -839,10 +826,10 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
avg_throat_diff += diff
count += 1
# Successfully connected this throat
poreTtbr.remove(throatIdx)
# Remove throat from neighbor
nborTtbr.remove(nborThroatIdx)
# Successfully connected this throat
ttbr.remove(throatIdx)
# remove fully-connected nbor from search neighborhood
if (len(nborTtbr) == 0):
......@@ -852,15 +839,12 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
cellList.expel(cpore)
# remove fully-connected pore from search neighborhood
if (len(ttbr) == 0):
if (len(poreTtbr) == 0):
cellList.expel(pore)
for cpore in copies[pore]:
cellList.expel(cpore)
# Free space
del cellList.poreMatchTable
# Finally, find next best candidates for left-out pores
# in REVERSED order
if (not mute): print("\nfinding next best candidates")
......@@ -872,21 +856,21 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
if (not mute):
print(f"progress {((it+1) / n)*100.:.1f}%", end="\r", flush=True)
ttbr = cellList.poreThroatIndices[poreIdx]
if (len(ttbr) == 0):
poreTtbr = ttbrs[poreIdx]
if (len(poreTtbr) == 0):
continue # Nothing left to do
# Expensive: Find alternative (multiprocessing)
candidates = cellList.find_candidates(poreIdx)
candidates = cellList.find_candidates(pore)
for throatIdx in ttbr.copy():
for throatIdx in poreTtbr.copy():
r = cellList.throatR[throatIdx]
targetLen = cellList.throatL[throatIdx]
r = Network.throatR[throatIdx]
targetLen = Network.throatL[throatIdx]
if (len(candidates) == 0):
# no candidates left; give up
k += len(ttbr) # update num unrealized throats
k += len(poreTtbr) # update num unrealized throats
break
nborIdx, lt = min(candidates.items(), key=lambda il: abs(il[1] - targetLen))
......@@ -911,13 +895,14 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
for cpore in copies[nbor]:
if cpore.index in candidates:
del candidates[cpore.index]
nborTtbr = ttbrs[nbor.index]
# connect pore to nbor
throats.append([pore, nbor, lbl, r])
nborTtbr = cellList.poreThroatIndices[nbor.index]
# remove most similar throat of nbor
ranking = [(nborThroatIdx, abs(cellList.throatL[nborThroatIdx] - lt)) \
ranking = [(nborThroatIdx, abs(Network.throatL[nborThroatIdx] - lt)) \
for nborThroatIdx in nborTtbr]
nborThroatIdx, _ = min(ranking, key=lambda il: il[1])
......@@ -926,7 +911,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
count += 1
# Remove throat from pore
ttbr.remove(throatIdx)
poreTtbr.remove(throatIdx)
# Remove throat from neighbor
nborTtbr.remove(nborThroatIdx)
......@@ -938,7 +923,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
cellList.expel(cpore)
# remove fully connected pore from search neighborhood
if (len(ttbr) == 0):
if (len(poreTtbr) == 0):
cellList.expel(pore)
for cpore in copies[pore]:
......@@ -985,7 +970,7 @@ 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, index=runningIndex, id=pore.id)
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
......@@ -996,7 +981,7 @@ def __add_buffer_layers(pores: List[Pore], targetsize: List[float],
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, index=runningIndex, id=pore.id)
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
......@@ -1202,7 +1187,7 @@ def save_network_to(filename: str, network: Network):
data=wrk)
def load_network_from(filename: str) -> Network:
def load_network_from(filename: str, isGenerator: bool = False) -> Network:
"""Load pore network from hdf5 file."""
import h5py
import numpy as np
......@@ -1229,13 +1214,17 @@ def load_network_from(filename: str) -> Network:
for id, strg in zip(f['pores/label/id'],f['pores/label/strg']):
pores[id].label = strg.decode('utf-8')
connector_method = network.connect_pores
if isGenerator:
connector_method = network.connect_pores_generator
# throats
tid = np.array(f['throats/id'])
throats = {} # id:throat dictionary
for id, pore1, pore2 in zip(tid,
[pores[id] for id in np.array(f['throats/pore1'])],
[pores[id] for id in np.array(f['throats/pore2'])]):
throats[id] = network.connect_pores(pore1, pore2, 0.0, id)
throats[id] = connector_method(pore1, pore2, 0.0, id)
# radius
for id, r in zip(tid, np.array(f['throats/r'])): throats[id].r = r
# labels
......
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