Commit da7acb9b authored by sfritschi's avatar sfritschi
Browse files

Further changes to generation algorithm

parent 43a159d8
......@@ -16,7 +16,7 @@ def main():
print("Min. number of throats per pore: %d" % min([len(pore.throats) for pore in basenet.pores]))
print("Avg. number of throats per pore: %f" % (sum([len(pore.throats) for pore in basenet.pores]) / len(basenet.pores)))
target = [2, 2, 1]
target = [2, 2, 2]
cutoff = 0.5 * max([basenet.ub[i] - basenet.lb[i] \
for i in range(len(basenet.ub))])
......
gen/plots/dendro.png

148 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
......@@ -100,6 +100,88 @@ class Network:
pore1.throats.add(throat); pore2.throats.add(throat)
return throat
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: Pore) -> List[int]:
from math import floor
# Shift pore position to be positive first (buffer layer)
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 (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
def distance(p1: List[float], p2: List[float]):
"""Distance between points p1 and p2."""
......@@ -370,89 +452,6 @@ def generate_imperial(basenet: Network, targetsize: List[float], \
for throat in throats: network.connect_pores(pore1=pores[throat[0]],
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: 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 (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
def generate_dendrogram(basenet: Network, targetsize: List[int], \
cutoff: float = float('inf'), sd: int = None, mute: bool = False) \
......@@ -548,46 +547,6 @@ 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)
# helper routines for connecting pores
# triangulation of pores excluding fully connected ones
def triangulation(pores):
# list of permissible pores (not fully connected ones)
p = [k for k, pore in enumerate(pores) if len(pore.throats) > 0]
neighborhood = {} # {pore:{neighboring pores}}
if (len(p) > 3):
# triangulation of permissible pores
import numpy as np
from scipy.spatial import Delaunay
tri = Delaunay(np.array([pores[k].pos for k in p])) # coplanar pnts!
# translate tri (point indices) into pore neighborhood dict
for j, k in enumerate(p):
nl = tri.vertex_neighbor_vertices[1][ \
tri.vertex_neighbor_vertices[0][j]: \
tri.vertex_neighbor_vertices[0][j+1]] # neighbor list
neighborhood[pores[k]] = {pores[p[n]] for n in nl}
elif (len(p) > 1):
for k in p:
neighborhood[pores[k]] = {pores[n] for n in p if n != k}
return neighborhood
# remove outsider from front, inside, and neighborhood
def expel(pore, outsider, neighborhood, front, inside, nnbrs):
# grow front by neighbors of outsider and shrink inside
if (outsider in front):
for nbor in neighborhood[outsider]:
if (nbor == pore) or (nbor in inside): continue # nbor is inside
l = distance(pore.pos, nbor.pos)
front[nbor] = l; inside[nbor] = l
del(front[outsider])
if (outsider in inside): del(inside[outsider])
# remove outsider from neighborhood
for nbor in neighborhood[outsider]:
nbors = neighborhood[nbor]
neighborhood[nbor] = nbors | (neighborhood[outsider] - {nbor})
neighborhood[nbor].remove(outsider)
nnbrs = max(nnbrs, len(neighborhood[nbor])) # update max. number of neighbors
del(neighborhood[outsider])
return nnbrs
"""
DEBUG: CELL LIST STATISTICS
......@@ -611,7 +570,6 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
if (not mute): print("\b"*21 + "connecting")
throats = []
k = 0 # count of unrealised throats
#nnbrs0 = 0; nnbrs = 0; # (initial) maximum number of neighbors
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)
......@@ -646,9 +604,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
if (not p in conpores)]
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
lt, nbor = min(ranking, key=lambda lp: abs(lp[0] - lt))
nborc = nbor # memorize potential periodic copy
# find original of buffer layer pore
if (nbor.label != LABELS[0]):
......@@ -660,12 +616,11 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
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])
_, nbor_throat = min(ranking, key=lambda lth: abs(lth[0] - lt))
nbor.throats.remove(nbor_throat)
# remove fully connected nbor from search neighborhood
if (len(nbor.throats) == 0):
# remove fully-connected nbor pore from candidates
......@@ -674,6 +629,9 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
poresSorted[nborIdx].remove(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)
# remove fully connected pore from search neighborhood
......
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