File wellmet/dicebox/circumtri.py changed (mode: 100644) (index f724535..5583667) |
... |
... |
class _CircumTri(_Exploration): |
426 |
426 |
|
|
427 |
427 |
def perform_topological_analysis(bx): |
def perform_topological_analysis(bx): |
428 |
428 |
assert bx.tri_space == 'G' |
assert bx.tri_space == 'G' |
429 |
|
bx.ta = sx.TopologyAnalysis(bx.convex_hull, bx.G, bx.Tri.tri.simplices, |
|
430 |
|
bx.Tri.tri.neighbors, bx.failsi) |
|
|
429 |
|
# tn_scheme, points, simplices, neighbors, failsi, facets, normals |
|
430 |
|
bx.ta = sx.TopologyAnalysis(bx.convex_hull.tn_scheme, bx.G, |
|
431 |
|
bx.Tri.tri.simplices, bx.Tri.tri.neighbors, |
|
432 |
|
bx.failsi, bx.convex_hull.simplices, |
|
433 |
|
bx.convex_hull.A) |
431 |
434 |
print("start integration of the convex envelope") |
print("start integration of the convex envelope") |
432 |
435 |
bx.ta.integrate_convex_envelope() |
bx.ta.integrate_convex_envelope() |
433 |
436 |
print("start integration of the internal walls") |
print("start integration of the internal walls") |
File wellmet/simplex.py changed (mode: 100644) (index a4bd4a4..1198138) |
... |
... |
from scipy import stats |
12 |
12 |
|
|
13 |
13 |
import quadpy |
import quadpy |
14 |
14 |
|
|
15 |
|
from collections import namedtuple, defaultdict |
|
|
15 |
|
from collections import namedtuple |
16 |
16 |
|
|
17 |
17 |
|
|
18 |
18 |
|
|
|
... |
... |
class CircumCenter: |
254 |
254 |
|
|
255 |
255 |
|
|
256 |
256 |
|
|
|
257 |
|
#č vlastní slovník se dycky hodí |
|
258 |
|
#č ten, na rozdil od defaultdict'a, neuklada chybejicí složky do slovníku |
|
259 |
|
class DefaultDict(dict): |
|
260 |
|
def __init__(self, nvar): |
|
261 |
|
self.nvar = nvar |
|
262 |
|
|
|
263 |
|
def __missing__(self, key): |
|
264 |
|
return np.zeros(self.nvar, dtype=float) |
|
265 |
|
|
|
266 |
|
|
257 |
267 |
|
|
258 |
268 |
|
|
259 |
269 |
# np.count_nonzero |
# np.count_nonzero |
|
... |
... |
class TopologyAnalysis: |
261 |
271 |
|
|
262 |
272 |
|
|
263 |
273 |
|
|
264 |
|
def __init__(self, convex_hull, points, simplices, neighbors, failsi): |
|
265 |
|
assert convex_hull.space == 'G' |
|
|
274 |
|
def __init__(self, tn_scheme, points, simplices, neighbors, failsi, facets, normals): |
266 |
275 |
self.simplices = simplices |
self.simplices = simplices |
267 |
276 |
self.neighbors = neighbors |
self.neighbors = neighbors |
268 |
277 |
self.failsi = failsi |
self.failsi = failsi |
269 |
|
self.convex_hull = convex_hull |
|
|
278 |
|
self.facets = facets |
|
279 |
|
self.normals = normals |
270 |
280 |
self.points = points |
self.points = points |
271 |
|
self.tn_scheme = convex_hull.tn_scheme |
|
|
281 |
|
self.tn_scheme = tn_scheme |
272 |
282 |
|
|
273 |
283 |
self.to_parse = set() |
self.to_parse = set() |
274 |
284 |
|
|
275 |
285 |
self.internal_walls = set() |
self.internal_walls = set() |
276 |
286 |
#self.external_facets = dict() |
#self.external_facets = dict() |
277 |
|
self.simplex_values = defaultdict(float) |
|
278 |
287 |
|
|
279 |
|
self.nvar = self.convex_hull.sample.nvar |
|
|
288 |
|
nsim, nvar = points.shape |
|
289 |
|
self.nvar = nvar |
|
290 |
|
self.nsim = nsim |
|
291 |
|
self.simplex_values = DefaultDict(nvar) |
|
292 |
|
|
|
293 |
|
|
|
294 |
|
|
280 |
295 |
|
|
281 |
296 |
|
|
282 |
297 |
def integrate_convex_envelope(self): |
def integrate_convex_envelope(self): |
283 |
|
facets = np.sort(self.convex_hull.simplices[:, ::-1], kind='mergesort', axis=1) |
|
|
298 |
|
facets = np.sort(self.facets[:, ::-1], kind='mergesort', axis=1) |
284 |
299 |
vertices_model = self.points[facets].transpose((1, 0, 2)) |
vertices_model = self.points[facets].transpose((1, 0, 2)) |
285 |
300 |
result = self.tn_scheme.integrate(self._norm_inside_callback, vertices_model) |
result = self.tn_scheme.integrate(self._norm_inside_callback, vertices_model) |
286 |
301 |
|
|
|
... |
... |
class TopologyAnalysis: |
291 |
306 |
|
|
292 |
307 |
#č nejdřív rouška - prointegrovalo se to vůbec kladně? |
#č nejdřív rouška - prointegrovalo se to vůbec kladně? |
293 |
308 |
#č pak normala - ona je pro tyhle ortogonálky konstantní |
#č pak normala - ona je pro tyhle ortogonálky konstantní |
|
309 |
|
result *= mask |
|
310 |
|
result *= self.normals.T |
294 |
311 |
|
|
295 |
|
normals = self.convex_hull.A |
|
296 |
312 |
|
|
297 |
|
result *= normals.T |
|
298 |
|
facet_integrals = np.average(result, weights=mask, axis=0) |
|
|
313 |
|
self.p_inside = np.sum(result, axis=1) |
299 |
314 |
|
|
300 |
|
self.p_inside = np.sum(facet_integrals) |
|
301 |
|
|
|
302 |
|
self.external_facets = dict(zip((tuple(facet) for facet in facets), facet_integrals)) |
|
|
315 |
|
self.external_facets = dict(zip((tuple(facet) for facet in facets), result.T)) |
303 |
316 |
|
|
304 |
317 |
|
|
305 |
318 |
|
|
|
... |
... |
class TopologyAnalysis: |
400 |
413 |
def _integrate_wall(self, vertices_model, inner_node): |
def _integrate_wall(self, vertices_model, inner_node): |
401 |
414 |
|
|
402 |
415 |
vectors = vertices_model[1:] - vertices_model[0] |
vectors = vertices_model[1:] - vertices_model[0] |
403 |
|
basis, __ = np.linalg.qr(vertices_model.T, mode='complete') |
|
|
416 |
|
basis, __ = np.linalg.qr(vectors.T, mode='complete') |
404 |
417 |
|
|
405 |
418 |
normal = basis[:, -1] |
normal = basis[:, -1] |
406 |
419 |
|
|
|
... |
... |
class TopologyAnalysis: |
440 |
453 |
#č nejdřív rouška - prointegrovalo se to vůbec kladně? |
#č nejdřív rouška - prointegrovalo se to vůbec kladně? |
441 |
454 |
#č pak normala - ona je pro tyhle ortogonálky konstantní |
#č pak normala - ona je pro tyhle ortogonálky konstantní |
442 |
455 |
|
|
|
456 |
|
|
|
457 |
|
values *= mask |
443 |
458 |
values *= normal |
values *= normal |
444 |
|
probability = vol * np.average(values, weights=mask) |
|
|
459 |
|
values *= vol * normal_sign |
445 |
460 |
|
|
446 |
|
return probability * normal_sign |
|
|
461 |
|
return values |
447 |
462 |
|
|
448 |
463 |
|
|
449 |
464 |
|
|
|
... |
... |
class TopologyAnalysis: |
492 |
507 |
|
|
493 |
508 |
|
|
494 |
509 |
event_id = self.event_ids[start_id] |
event_id = self.event_ids[start_id] |
|
510 |
|
|
|
511 |
|
nvar = self.nvar |
495 |
512 |
|
|
496 |
|
region = Region() |
|
|
513 |
|
region = Region(nvar) |
497 |
514 |
region.event_id = event_id |
region.event_id = event_id |
498 |
515 |
|
|
499 |
516 |
# 0=success, 1=failure, 2=mix |
# 0=success, 1=failure, 2=mix |
|
... |
... |
class TopologyAnalysis: |
505 |
522 |
self.greens.append(region) |
self.greens.append(region) |
506 |
523 |
|
|
507 |
524 |
|
|
508 |
|
nvar = self.nvar |
|
|
525 |
|
|
509 |
526 |
wall_indices = np.empty(nvar, dtype=int) |
wall_indices = np.empty(nvar, dtype=int) |
510 |
527 |
vertex_arange = np.arange(nvar + 1) |
vertex_arange = np.arange(nvar + 1) |
511 |
528 |
|
|
|
... |
... |
class TopologyAnalysis: |
547 |
564 |
class Region: |
class Region: |
548 |
565 |
__slots__ = ('simplex_index', 'neighbors', 'probability', 'event_id') |
__slots__ = ('simplex_index', 'neighbors', 'probability', 'event_id') |
549 |
566 |
|
|
550 |
|
def __init__(self): |
|
|
567 |
|
def __init__(self, nvar): |
551 |
568 |
self.simplex_index = set() |
self.simplex_index = set() |
552 |
569 |
self.neighbors = set() |
self.neighbors = set() |
553 |
|
self.probability = 0 |
|
|
570 |
|
self.probability = np.zeros(nvar) |
554 |
571 |
|
|
555 |
572 |
# def __repr__(self): |
# def __repr__(self): |
556 |
573 |
# return "%s(**%s)"%(self.__class__.__name__, repr(self.init_parameters())) |
# return "%s(**%s)"%(self.__class__.__name__, repr(self.init_parameters())) |