Commit 7d243bad authored by sfritschi's avatar sfritschi
Browse files

Fixed issue with periodic copies if targetsize contains 1

parent 1c2cc8af
......@@ -20,7 +20,7 @@ def main():
print("Min. number of throats per pore: %d" % min(pore_throat_counts))
print("Avg. number of throats per pore: %f" % (sum(pore_throat_counts) / len(pore_throat_counts)))
target = [2, 2, 2]
target = [5, 5, 5]
print("Target size: {}".format(target))
cutoff = 0.5 * max([basenet.ub[i] - basenet.lb[i] \
for i in range(len(basenet.ub))])
......
......@@ -107,8 +107,9 @@ class CellList:
find_candidates() assumes pore is located in interior of domain (non-periodic)
"""
def __init__(self, pores: List[Pore], domainSize: List[float], \
cellSize: float, basenet: Network, nInterior: int):
def __init__(self, pores: List[Pore], copies: Dict[Pore, Set[Pore]], \
domainSize: List[float], cellSize: float, \
basenet: Network, nInterior: int, isPeriodicCheck: bool):
from math import ceil
self.dim = len(domainSize)
self.Lmax = basenet.Lmax
......@@ -128,9 +129,18 @@ class CellList:
# Offset in pore-throat table for given pore index
self.throatOffsets = [0] # size: #pores + 1 (non-periodic)
offset = 0
# Used to keep track of realized throats
# Used to keep track of already realized throats
self.poreThroatIndices = [set() for _ in range(nInterior)]
# Additional information needed for correctness if targetsize contains 1
self.isPeriodicCheck = isPeriodicCheck
if isPeriodicCheck:
# Maps periodic copy index to respective index inside domain
# {int: int}
self.originalIdx = {}
# Maps index from inside domain to List of all its periodic copies
# {int: List[int]}
self.periodicCopies = {}
L = [ub-lb for lb,ub in zip(basenet.lb,basenet.ub)]
# Sort pores according to cell membership
self.poresSorted = [set() for _ in range(self.totalCells)]
......@@ -146,14 +156,28 @@ class CellList:
offset += nthroats
self.throatOffsets.append(offset)
if isPeriodicCheck:
self.originalIdx[poreIdx] = poreIdx # identity
copyIndices = []
for cpore in copies[pore]:
copyIndices.append(cpore.index)
self.periodicCopies[poreIdx] = copyIndices
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
elif isPeriodicCheck: # periodic copy of pore in interior
# Determine index of original pore
originalIdx, _ = pore.label.split(' ',1)
self.originalIdx[poreIdx] = int(originalIdx)
self.poreThroatTable = mp.RawArray('i', throatIdx)
throatTotal = throatIdx
self.poreThroatTable = mp.RawArray('i', throatTotal)
# Helper for converting 3D index triple into 1D (flattened) index
def flatten(self, i: int, j: int, k: int) -> int:
......@@ -198,9 +222,19 @@ class CellList:
matchIdx, matchLen = min(candidates.items(), key=lambda il: abs(il[1] - lt))
# Remove match from candidates
# If targetsize contains 1 need to remove periodic copies as well
del candidates[matchIdx]
# Remove periodic copies as well if targetsize contains 1
if self.isPeriodicCheck:
# Get index of original pore in interior of domain
originalIdx = self.originalIdx[matchIdx]
if originalIdx in candidates:
del candidates[originalIdx]
for cporeIdx in self.periodicCopies[originalIdx]:
if cporeIdx in candidates:
del candidates[cporeIdx]
self.set_throat(poreIdx, matchIdx, idx)
# Compute dictionary of all nbor candidates of given pore
......@@ -685,7 +719,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
"""
import multiprocessing as mp
from random import seed, random, randint
from itertools import product
from itertools import product, accumulate
seed(sd)
d = len(targetsize) # number of spatial dimensions
# check target size multiplicator (must be >= 1)
......@@ -772,11 +806,12 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
# Flag for handling the case where original nbor pore and its periodic
# copies are potential candidates of a given pore at the same time
# (necessary only if target size is 1 in at least one dimension)
#isPeriodicCheck = 1 in targetsize
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, trueDomainSize, basenet.Lmax, basenet, n)
cellList = CellList(pores, copies, trueDomainSize, basenet.Lmax, \
basenet, n, isPeriodicCheck)
# Evenly distribute pores among threads
remainder = n % nthreads
......@@ -784,14 +819,13 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
for i in range(remainder):
loads[i] += 1
from itertools import accumulate
# Compute inclusive-scan
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 = range(displs[tid], displs[tid+1])
worker = mp.Process(target=cellList.match_maker, args=(poreList, tid))
......@@ -827,10 +861,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
if (throatIdx not in ttbr):
continue # Throat already realized
r = cellList.throatR[throatIdx]
targetLen = cellList.throatL[throatIdx]
nbor = pores[match]
lt = distance(pore.pos, nbor.pos)
nborc = nbor # memorize potential periodic copy
# find original of buffer layer pore
......@@ -841,12 +872,17 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
lbl = LABELS[0]
nborTtbr = cellList.poreThroatIndices[nbor.index]
total += 1
if (len(nborTtbr) == 0):
n_already_taken += 1
continue
r = cellList.throatR[throatIdx]
targetLen = cellList.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
......@@ -863,14 +899,14 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
# Successfully connected this throat
ttbr.remove(throatIdx)
# remove fully connected nbor from search neighborhood
# remove fully-connected nbor from search neighborhood
if (len(nborTtbr) == 0):
cellList.expel(nbor)
for cpore in copies[nbor]:
cellList.expel(cpore)
# remove fully connected pore from search neighborhood
# remove fully-connected pore from search neighborhood
if (len(ttbr) == 0):
cellList.expel(pore)
......@@ -904,9 +940,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
break
nborIdx, lt = min(candidates.items(), key=lambda il: abs(il[1] - targetLen))
del candidates[nborIdx]
nbor = pores[nborIdx]
nborc = nbor # memorize potential periodic copy
# find original of buffer layer pore
......@@ -916,6 +950,18 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
else:
lbl = LABELS[0]
# 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]
# connect pore to nbor
throats.append([pore, nbor, lbl, r])
......
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