iam-git / WellMet (public) (License: MIT) (since 2021-08-31) (hash sha1)
WellMet is pure Python framework for spatial structural reliability analysis. Or, more specifically, for "failure probability estimation and detection of failure surfaces by adaptive sequential decomposition of the design domain".
Clone URLs: https://rocketgit.com/user/iam-git/WellMet ssh://rocketgit@ssh.rocketgit.com/user/iam-git/WellMet git://git.rocketgit.com/user/iam-git/WellMet
release trunk
List of commits:
Subject Hash Author Date (UTC)
simplex: make _Triangulation more clean; rework CubatureIntegration fb7dd9fef8592944fbd93a3f07b44c2355d45e20 I am 2022-12-28 15:48:20
introduce shell module, clean version of ghull b13e9b2b15d109770e0c3fbdd5bb7b2a3b5bf5fc I am 2022-12-28 15:45:27
qt_gui.qt_plot.CandidatesWidget: regression fix a75874b5a53e7b420672b582f3ff415b17a32e9b I am 2022-12-23 12:19:51
qt_gui.qt_pairwise.ContactWidget: add mask function 6f88fec7db92e5cd21826a06f42a609a3e03cb6b I am 2022-12-21 15:20:01
qt_gui.qt_pairwise.ContactWidget: add option for mixed pairs only adjacency search. 2f87e8a18908de607b8abc0779448ca3095c134c I am 2022-12-21 14:30:23
qt_gui.qt_dicebox: reflect changes in CircumTri class 9821d389240a383f0b13b7dd92cd501e1c40c537 I am 2022-12-15 09:36:16
dicebox.circumtri: employ _Exploration class. Keep naive CircumTri class in renamed module b8605fb24d55b3296ceac34afba541a5b7ee3070 I am 2022-12-15 09:33:59
convex_hull.QHull.get_exploration_vector: self -> hull fix dae47ef37b5871a98d2ea274cc35aeb2a4eebe1c I am 2022-12-15 09:30:31
qt_gui/qt_dicebox: add DumpExplorationWidget b687895f06c780418a12eea59bd052ad7cd348f3 I am 2022-12-15 01:31:14
dicebox: add helper _exploration module 836ae6cd12f1a04f1e2fc3423e46777d16ebad7d I am 2022-12-15 01:30:20
convex_hull: add get_exploration_vector() method 795259d91b64c732e60a83a1d9e3b89c2141b003 I am 2022-12-15 01:28:12
qt_gui.qt_gui: dot not create graph widgets if estimations are not present f070ed140a135a980e4436ed3779455c0e1bf399 I am 2022-12-14 02:25:45
qt_gui.gl_plot.CandidatesWidget: apply changes for GLplot too 03550040b60c67d247e26425e35fac565dd88886 Aleksei Gerasimov 2022-12-14 01:56:01
qt_gui.qt_plot.CandidatesWidget: add CandyNodes support 2357cdd92170b1f660b15c92aec695d3f7c108b3 I am 2022-12-14 01:08:30
qt_gui.qt_graph_widgets.EstimationGraph: reimplement piece of zerosafe logic 342a1b7bfddb7ba5695b3f0e017c7f75811db1f0 I am 2022-12-13 22:15:54
qt_gui.qt_graph_widgets: update x range when box runned c52052a9c3f6926d03e8c39f7642aa0a8864fdc3 I am 2022-12-13 10:23:45
qt_gui.qt_graph_widgets.EstimationGraph: add outside and mixed lines df6c06e4287492855d613b1ae46edfca6a4757d3 I am 2022-12-12 06:13:55
qt_gui: replace old copypaste by cleaned up EstimationGraph 515d0d8a892984f8a7667da874972daff3318d74 I am 2022-12-12 05:36:43
qt_gui: prepare new BoxEstimationData class; rework ErrorGraph 37387100715f980771a8a37f99eec3e723759cd3 I am 2022-12-11 23:17:06
dicebox.circumtri: replace ndim by nvar to make TriEstimation pandas-friendly 0824222ae4c9c1f3fede971481816ed9ded9707a I am 2022-12-11 23:14:05
Commit fb7dd9fef8592944fbd93a3f07b44c2355d45e20 - simplex: make _Triangulation more clean; rework CubatureIntegration
Author: I am
Author date (UTC): 2022-12-28 15:48
Committer name: I am
Committer date (UTC): 2022-12-28 15:48
Parent(s): b13e9b2b15d109770e0c3fbdd5bb7b2a3b5bf5fc
Signer:
Signing key:
Signing status: N
Tree: 1a84c5ffbbb8720360e595a61769bfe005c96fa3
File Lines added Lines deleted
wellmet/simplex.py 279 179
File wellmet/simplex.py changed (mode: 100644) (index 945599e..97a9d84)
... ... from scipy import stats
12 12
13 13 import quadpy import quadpy
14 14
15 from collections import namedtuple
16
17
18
15 19 #č napadlo mě zababáchnuť třidu, která by se sama starala o všem co se tyče #č napadlo mě zababáchnuť třidu, která by se sama starala o všem co se tyče
16 20 #č vnější domény. Nešlo mě totíž to udělat jednou funkcí, bylo by velmi #č vnější domény. Nešlo mě totíž to udělat jednou funkcí, bylo by velmi
17 21 #č špatné z hlediska zodpovednosti kódu. Tak to všecko zabalíme to třidy #č špatné z hlediska zodpovednosti kódu. Tak to všecko zabalíme to třidy
 
... ... class CircumCenter:
246 250 x = np.linalg.solve(self.A, self.b) x = np.linalg.solve(self.A, self.b)
247 251 return vertices.T @ x[:-1] return vertices.T @ x[:-1]
248 252
253
249 254 #č šablona třidy #č šablona třidy
250 255 class _Triangulation: class _Triangulation:
251 def setup(sx, sample_box, tri_space='Rn', issi=None, weighting_space=None, \
252 incremental=True, on_add_simplex=None, on_delete_simplex=None):
256 def tri_setup(sx, sample_box, tri_space='Rn', incremental=True):
253 257
254 258 sx.sample_box = sample_box sx.sample_box = sample_box
255 259 sx.tri_space = tri_space sx.tri_space = tri_space
256 sx.PDF = sample_box.pdf(tri_space)
257 sx.failsi = sample_box.failsi
258 sx.issi = issi #č ISSI potřebujem pro tri_estimation
259
260 if weighting_space is None:
261 sx.weighting_space = tri_space
262 else:
263 sx.weighting_space = weighting_space
264 260
265 #č kolbeky
266 sx.on_add_simplex = on_add_simplex
267 sx.on_delete_simplex = on_delete_simplex
268
269 #оӵ кылсузъет кылдытом
270 sx.simplex_stats = dict()
271 #č ty množiny jsou tak trošku overkill, ale budiž
272 261 sx.simplices_set = set() sx.simplices_set = set()
273 262
274 263 sx.newly_estimated = 0 sx.newly_estimated = 0
 
... ... class _Triangulation:
278 267 #č tri - Deloneho triangulace #č tri - Deloneho triangulace
279 268 #č žádné chyby nechytám #č žádné chyby nechytám
280 269 #čs když se tringulace nepovede tak nemáme čo robiť #čs když se tringulace nepovede tak nemáme čo robiť
281 # incremental triangulation require one more point
270 # incremental triangulation requires one more point
282 271 tri_plan = getattr(sample_box, tri_space) tri_plan = getattr(sample_box, tri_space)
283 272 sx.tri = spatial.Delaunay(tri_plan, incremental=incremental) sx.tri = spatial.Delaunay(tri_plan, incremental=incremental)
284 273 if len(sx.tri.coplanar): if len(sx.tri.coplanar):
 
... ... class _Triangulation:
305 294 return np.int8(np.where(has_failure, np.where(all_failure, 1, 2), 0)) return np.int8(np.where(has_failure, np.where(all_failure, 1, 2), 0))
306 295
307 296
297 def get_nfailures(sx, simplices=None):
298 if simplices is None:
299 simplices = sx.tri.simplices
300
301 in_failure = np.isin(simplices, sx.sample_box.failure_points)
302 return np.sum(in_failure, axis=1)
303
304
308 305
309 306 #č tato funkce běží 91% času #č tato funkce běží 91% času
310 307 # bottleneck function # bottleneck function
311 308 def update(sx): def update(sx):
309 simplices_set_to_estimate = sx._update()
310 gc.collect()
311
312 sx.estimate_simplices(simplices_set_to_estimate)
313
314
315 def _update(sx):
312 316 """ """
313 317 Triangulace zajistěně existuje Triangulace zajistěně existuje
314 318 """ """
 
... ... class _Triangulation:
322 326 #č reference, jenom dejte vědet, #č reference, jenom dejte vědet,
323 327 #č že máme triangulaci aktualizovat #č že máme triangulaci aktualizovat
324 328
325 former_simplices = sx.tri.simplices
326 former_nsimplex = sx.tri.nsimplex
329 sx._tri_update() #č jako vždy, chyby nechytáme
327 330
328 tri_plan = getattr(sx.sample_box, sx.tri_space)
329 #č jako vždy, chyby nechytáme
330 #sx.tri.add_points(tri_plan[sx.tri.npoints:])
331 sx._tri_update(tri_plan)
332
333 if len(sx.tri.coplanar): # pokud triangulace není v pořadku
334 #print('triangulace v pořádku není')
335 print('Triangulation has coplanar points:', sx.tri.coplanar)
336
337 331
338 332 #č vyhodit ty pomalé pytloviny, co tu byly #č vyhodit ty pomalé pytloviny, co tu byly
339 333 #č (tady bylo něco jako tringulace - 1,45 s, drbání kolem - 1366 s), #č (tady bylo něco jako tringulace - 1,45 s, drbání kolem - 1366 s),
 
... ... class _Triangulation:
372 366 #č necháme zbytek jednotlivým podtřídám #č necháme zbytek jednotlivým podtřídám
373 367 #č co jsem viděl, voláme tyhle funkce jenom my #č co jsem viděl, voláme tyhle funkce jenom my
374 368 sx._invalidate_simplices(to_invalidate) sx._invalidate_simplices(to_invalidate)
375 gc.collect()
376 sx._estimate_simplices(to_estimate)
377
378 369
379 def _tri_update(sx, tri_plan):
380 #č separujeme, abychom vědeli, kolik času žere samotný QHull
381 sx.tri.add_points(tri_plan[sx.tri.npoints:])
382 #č ale jo, nejen QHull
383 sx.PDF = sx.sample_box.pdf(sx.tri_space)
384 sx.failsi = sx.sample_box.failsi
385 370
386
387
388 def _estimate_simplices(sx, simplices_set_to_estimate):
389 371 #č zde spolehám na to, že pořadí indexů se nikdy nezmění #č zde spolehám na to, že pořadí indexů se nikdy nezmění
390 372 #č tj. [12, 13, 26] najednou nestane [26, 12, 13] #č tj. [12, 13, 26] najednou nestane [26, 12, 13]
391 373 #č (dá se něco takovýho očekavát podle toho co jsem čet v dokumentaci) #č (dá se něco takovýho očekavát podle toho co jsem čet v dokumentaci)
392
374 #
393 375 # here "simplices_set_to_estimate" is a set of tuples # here "simplices_set_to_estimate" is a set of tuples
394 simplices = np.array(list(simplices_set_to_estimate))
395 for simplex in simplices:
396 sx.estimate_simplex(simplex)
397
376 #simplices_to_estimate = np.array(list(simplices_set_to_estimate))
377 return to_estimate
378
379
380 def _tri_update(sx):
381 tri_plan = getattr(sx.sample_box, sx.tri_space)
382 #č separujeme, abychom vědeli, kolik času žere samotný QHull
383 sx.tri.add_points(tri_plan[sx.tri.npoints:])
384
385 if len(sx.tri.coplanar): # pokud triangulace není v pořadku
386 #print('triangulace v pořádku není')
387 print('Triangulation has coplanar points:', sx.tri.coplanar)
398 388
399 389
400 390 def is_mixed(bx, simplices=None): def is_mixed(bx, simplices=None):
 
... ... class _FullTriangulation(_Triangulation):
449 439 sx.on_delete_simplex(indices=np.array(simplex)) sx.on_delete_simplex(indices=np.array(simplex))
450 440
451 441
442 def estimate_simplices(sx, simplices_set_to_estimate):
443 for simplex in simplices_set_to_estimate:
444 sx.estimate_simplex(np.array(simplex))
445
446
452 447 #č vyhodil jsem simplex_id'y #č vyhodil jsem simplex_id'y
453 448 def estimate_simplex(sx, indices): def estimate_simplex(sx, indices):
454 449 # -1 = 'outside', 0=success, 1=failure, 2=mix # -1 = 'outside', 0=success, 1=failure, 2=mix
 
... ... class _FastTriangulation(_Triangulation):
546 541 sx.on_delete_simplex(indices=np.array(simplex)) sx.on_delete_simplex(indices=np.array(simplex))
547 542
548 543
544 def estimate_simplices(sx, simplices_set_to_estimate):
545 for simplex in simplices_set_to_estimate:
546 sx.estimate_simplex(np.array(simplex))
547
548
549 549 #č vyhodil jsem simplex_id'y #č vyhodil jsem simplex_id'y
550 550 def estimate_simplex(sx, indices): def estimate_simplex(sx, indices):
551 551
 
... ... class _SamplingTriangulation:
577 577
578 578 sx.simplex_budget = simplex_budget sx.simplex_budget = simplex_budget
579 579 sx.design = design sx.design = design
580
581 sx.setup(sample_box, tri_space=tri_space, issi=issi,\
582 weighting_space=weighting_space, incremental=incremental,\
583 on_add_simplex=on_add_simplex,\
580
581 sx.tri_setup(sample_box, tri_space=tri_space, incremental=incremental)
582
583 sx.setup(issi=issi, weighting_space=weighting_space,
584 on_add_simplex=on_add_simplex,
584 585 on_delete_simplex=on_delete_simplex) on_delete_simplex=on_delete_simplex)
585 586
586 587
587 588
589 def setup(sx, issi=None, weighting_space=None,
590 on_add_simplex=None, on_delete_simplex=None):
591
592 sx.issi = issi #č ISSI potřebujem pro tri_estimation
593
594 if weighting_space is None:
595 sx.weighting_space = tri_space
596 else:
597 sx.weighting_space = weighting_space
598
599 #č kolbeky
600 sx.on_add_simplex = on_add_simplex
601 sx.on_delete_simplex = on_delete_simplex
602
603 #оӵ кылсузъет кылдытом
604 sx.simplex_stats = dict()
605
606
607
588 608 class _CubatureTriangulation: class _CubatureTriangulation:
589 609 def __init__(sx, sample_box, tn_scheme=None, tri_space='Rn', issi=None, weighting_space=None, \ def __init__(sx, sample_box, tn_scheme=None, tri_space='Rn', issi=None, weighting_space=None, \
590 610 incremental=True, on_add_simplex=None, on_delete_simplex=None): incremental=True, on_add_simplex=None, on_delete_simplex=None):
591 611
592 sx.setup(sample_box, tri_space=tri_space, issi=issi,\
593 weighting_space=weighting_space, incremental=incremental,\
594 on_add_simplex=on_add_simplex,\
612 sx.tri_setup(sample_box, tri_space=tri_space, incremental=incremental)
613 _SamplingTriangulation.setup(sx, issi=issi,
614 weighting_space=weighting_space,
615 on_add_simplex=on_add_simplex,
595 616 on_delete_simplex=on_delete_simplex) on_delete_simplex=on_delete_simplex)
596 617
597 618 #č Hodil by se nám levný odhad ve chvili, #č Hodil by se nám levný odhad ve chvili,
 
... ... class JustCubatureTriangulation(_FastTriangulation, _CubatureTriangulation):
1015 1036
1016 1037
1017 1038
1039 FastCubatureEstimation = namedtuple('FastCubatureEstimation', (
1040 'nvar',
1041 'nsim',
1042 'failure',
1043 'mix',
1044 'vertex_ratio_estimation',
1045 'vertex_estimation',
1046 'weighted_ratio_estimation',
1047 'nsimplex',
1048 #'tn_scheme',
1049 'tn_scheme_points',
1050 'newly_invalidated',
1051 'newly_estimated',
1052 'failure_simplices',
1053 'mixed_simplices',
1054 'ncoplanar'
1055 ))
1018 1056
1019 1057 #č zadavame v každem integračním bodě očekavaní pravděpodobnosti poruchy #č zadavame v každem integračním bodě očekavaní pravděpodobnosti poruchy
1020 1058 # p_f + p_s = 1 # p_f + p_s = 1
 
... ... class JustCubatureTriangulation(_FastTriangulation, _CubatureTriangulation):
1031 1069 # P_simplex = sum(f * w) # P_simplex = sum(f * w)
1032 1070 # Pf_simplex = pf_simplex * P_simplex = sum(p_f * w) * sum(f * w) # Pf_simplex = pf_simplex * P_simplex = sum(p_f * w) * sum(f * w)
1033 1071 # (p_f1 * w1 + p_f2*w2 + p_f3*w3) * (f1 * w1 + f2*w2 + f3*w3) # (p_f1 * w1 + p_f2*w2 + p_f3*w3) * (f1 * w1 + f2*w2 + f3*w3)
1034 class BetterCubatureIntegration(_Triangulation):
1035 def __init__(sx, sample_box, tn_scheme, tri_space='Rn', \
1036 incremental=True, on_add_simplex=None, on_delete_simplex=None):
1072 class CubatureIntegration(_Triangulation):
1073 def __init__(sx, sample_box, tn_scheme, tri_space='Rn', incremental=True,
1074 on_failure_added=lambda *__: None, on_mixed_added=lambda *__: None,
1075 on_delete_simplex=lambda __: None):
1037 1076
1038 sx.setup(sample_box, tri_space=tri_space, issi=None,\
1039 weighting_space=None, incremental=incremental,\
1040 on_add_simplex=on_add_simplex,\
1041 on_delete_simplex=on_delete_simplex)
1077 sx.tri_setup(sample_box, tri_space=tri_space, incremental=incremental)
1042 1078
1043 1079 sx.tn_scheme = tn_scheme sx.tn_scheme = tn_scheme
1044 1080
 
... ... class BetterCubatureIntegration(_Triangulation):
1046 1082 #č ptat se na něco u skříňky je extrémně dráho #č ptat se na něco u skříňky je extrémně dráho
1047 1083 sx.f = sample_box.f_model sx.f = sample_box.f_model
1048 1084
1085 sx.PDF = sample_box.pdf(tri_space)
1086 sx.failsi = sample_box.failsi
1087
1088 #č kolbeky
1089 sx.on_failure_added = on_failure_added
1090 sx.on_mixed_added = on_mixed_added
1091 sx.on_delete_simplex = on_delete_simplex
1092
1093 #оӵ кылсузъет кылдытом
1094 sx.failure_simplices = dict()
1095 sx.mixed_simplices = dict()
1096
1049 1097
1050 1098 def integrate(sx): def integrate(sx):
1051 #č Metoda musí simplexům přiřazovat jev
1052 # 0=success, 1=failure, 2=mix
1053 #č vyhodil jsem simplex_id'y
1054 event_ids = sx.get_events()
1099 simplices = sx.tri.simplices
1100 in_failure = np.isin(simplices, sx.sample_box.failure_points)
1055 1101
1056 # 0=success, 1=failure, 2=mix
1057 for simplex in sx.tri.simplices[event_ids == 1]:
1058 #č ty množiny jsou super
1059 sx.simplices_set.add(tuple(simplex))
1060 sx._integrate_red_simplex(simplex)
1061 for simplex in sx.tri.simplices[event_ids == 2]:
1102 has_failure = in_failure.any(axis=1)
1103 simplices = simplices[has_failure]
1104 in_failure = in_failure[has_failure]
1105
1106 all_failure = in_failure.all(axis=1)
1107 for indices in simplices[all_failure]:
1108 simplex = tuple(indices)
1062 1109 #č ty množiny jsou super #č ty množiny jsou super
1063 sx.simplices_set.add(tuple(simplex))
1064 sx._integrate_mixed_simplex(simplex)
1110 sx.simplices_set.add(simplex)
1111 sx._integrate_red_simplex(simplex, indices)
1065 1112
1113 simplices = simplices[~all_failure]
1114 in_failure = in_failure[~all_failure]
1115
1116 frs = np.sum(in_failure, axis=1) / (sx.sample_box.nvar + 1)
1117 pdfs = sx.PDF[simplices]
1118 fpdfs = np.sum(pdfs * in_failure, axis=1)
1119 wfrs = fpdfs / np.sum(pdfs, axis=1)
1120 mfpdfs = fpdfs / (sx.sample_box.nvar + 1)
1121
1122 for indices, fr, wfr, mfpdf in zip(simplices, frs, wfrs, mfpdfs):
1123 #č ty množiny jsou super
1124 simplex = tuple(indices)
1125 #č ty množiny jsou super
1126 sx.simplices_set.add(simplex)
1127 sx._integrate_mixed_simplex(simplex, indices, fr, wfr, mfpdf)
1066 1128
1129
1067 1130
1068 1131
1069 1132 #č vyhodil jsem simplex_id'y #č vyhodil jsem simplex_id'y
 
... ... class BetterCubatureIntegration(_Triangulation):
1073 1136
1074 1137 #č ty simplexy NEmusí tam být, #č ty simplexy NEmusí tam být,
1075 1138 for simplex in simplices_set_to_delete: for simplex in simplices_set_to_delete:
1076 sx.simplex_stats.pop(simplex, None)
1139 sx.failure_simplices.pop(simplex, None)
1140 sx.mixed_simplices.pop(simplex, None)
1077 1141
1078 if sx.on_delete_simplex is not None:
1079 #č Bacha, teď posíláme tuple.
1080 #č Kód na té straně úplně stejně bude vyhazovat
1081 #č záznamy ze slovníků
1082 sx.on_delete_simplex(simplex)
1142 #č Bacha, teď posíláme tuple.
1143 #č Kód na té straně úplně stejně bude vyhazovat
1144 #č záznamy ze slovníků
1145 sx.on_delete_simplex(simplex)
1083 1146
1084 1147
1085 #č vyhodil jsem simplex_id'y
1086 def estimate_simplex(sx, indices):
1148 def estimate_simplices(sx, simplices_set_to_estimate):
1149 #č aktualizace krajské urovně
1150 sx.PDF = PDF = sx.sample_box.pdf(sx.tri_space)
1151 sx.failsi = failsi = sx.sample_box.failsi
1087 1152
1088 #č zkusím funkci návrhnout tak, že
1089 1153 #ё вызывающая функция запускает estimate_simplex #ё вызывающая функция запускает estimate_simplex
1090 1154 #ё на всём подряд без разбору. #ё на всём подряд без разбору.
1091 1155 #č Našim úkolem je zjistit co je simplex zač #č Našim úkolem je zjistit co je simplex zač
1092 1156 #č a ty zelené ignorovat #č a ty zelené ignorovat
1093
1094 #č žádný intersect, žádný setdiff
1095 #č tohle je nejrychlejší!
1096 indices_failsi = sx.sample_box.failsi[indices]
1097 1157
1098 # fp like a failure points. Number of failure points
1099 fp = len(indices_failsi[indices_failsi]) # the fastest solution
1158 # fp is like failure points. Number of failure points
1159 max_fp = sx.sample_box.nvar + 1
1100 1160
1101 # -1 = 'outside', 0=success, 1=failure, 2=mix
1102 #č simplex -1 mít nemůže
1103 if fp == len(indices):
1104 sx._integrate_red_simplex(indices)
1105 elif fp > 0:
1106 sx._integrate_mixed_simplex(indices)
1161 for simplex in simplices_set_to_estimate:
1162 indices = np.array(simplex)
1163
1164 #č žádný intersect, žádný setdiff
1165 #č tohle je nejrychlejší!
1166 indices_failsi = failsi[indices]
1167
1168 # fp is like failure points. Number of failure points
1169 fp = len(indices_failsi[indices_failsi]) # the fastest solution
1170
1171 if fp == max_fp:
1172 sx._integrate_red_simplex(simplex, indices)
1173 elif fp > 0:
1174 fr = fp / max_fp
1175 pdf = PDF[indices]
1176 fpdf = np.sum(pdf[indices_failsi])
1177 # same as np.average(failsi, weights=pdf), but faster
1178 wfr = fpdf / np.sum(pdf)
1179 mfpdf = fpdf / max_fp
1180 sx._integrate_mixed_simplex(simplex, indices, fr, wfr, mfpdf)
1107 1181
1108 1182
1109 1183
1110 def _integrate_mixed_simplex(sx, indices):
1184 # def _integrate_mixed_simplex(sx, simplex, indices):
1185 #
1186 # vertices_model = sx.tri.points[indices]
1187 # PDF = sx.PDF[indices]
1188 # failsi = sx.failsi[indices]
1189 #
1190 # # quadpy
1191 # def _get_pdf(x):
1192 # dmv = 1 / spatial.distance.cdist(x.T, vertices_model)
1193 # dmw = dmv * PDF
1194 #
1195 # # same as np.sum(dmv[:, failsi], axis=1) / np.sum(dmv, axis=1)
1196 # pfv = np.sum(dmv.T[failsi], axis=0) / np.sum(dmv, axis=1)
1197 # pfw = np.sum(dmw.T[failsi], axis=0) / np.sum(dmw, axis=1)
1198 #
1199 # # side effect
1200 # if sx.on_add_simplex is not None:
1201 # nodes = sx.f.new_sample(x.T, sx.tri_space)
1202 # fx = nodes.pdf(sx.tri_space)
1203 # # -1=outside, 0=success, 1=failure, 2=mix
1204 # args = {'event_id':2, 'indices':indices}
1205 # sx.on_add_simplex(CandyNodes(nodes, args, pfv=pfv, pfw=pfw))
1206 # else:
1207 # fx = sx.f.sample_pdf(x.T, sx.tri_space)
1208 #
1209 # return fx, fx*pfv, fx*pfw
1210 #
1211 #
1212 # simplex_measures = sx.tn_scheme.integrate(_get_pdf, vertices_model)
1213 # #č místo fallback integrací
1214 # mask = simplex_measures < 0
1215 # if np.any(mask):
1216 # print("Negative measures have occured in simplex %s" % indices)
1217 # print("Rozbíla se nám integrace, integráční schema je na houby.")
1218 # simplex_measures[mask] = 0
1219 #
1220 # p_mixed, pfv_simplex, pfw_simplex = simplex_measures
1221 #
1222 #
1223 # assert p_mixed >= 0
1224 # assert pfv_simplex >= 0
1225 # assert pfw_simplex >= 0
1226 #
1227 # #č přídavám jednou, sčítám pořad dokola
1228 # simplex_estimation = np.array((p_mixed, 0, pfv_simplex, pfw_simplex))
1229 # sx.simplex_stats[simplex] = simplex_estimation
1230
1231
1232 def _integrate_mixed_simplex(sx, simplex, indices, fr, wfr, mfpdf):
1111 1233
1112 1234 vertices_model = sx.tri.points[indices] vertices_model = sx.tri.points[indices]
1113 PDF = sx.PDF[indices]
1114 failsi = sx.failsi[indices]
1115 1235
1116 # quadpy
1117 def _get_pdf(x):
1118 dmv = 1 / spatial.distance.cdist(x.T, vertices_model)
1119 dmw = dmv * PDF
1120
1121 # same as np.sum(dmv[:, failsi], axis=1) / np.sum(dmv, axis=1)
1122 pfv = np.sum(dmv.T[failsi], axis=0) / np.sum(dmv, axis=1)
1123 pfw = np.sum(dmw.T[failsi], axis=0) / np.sum(dmw, axis=1)
1124
1125 # side effect
1126 if sx.on_add_simplex is not None:
1127 nodes = sx.f.new_sample(x.T, sx.tri_space)
1128 fx = nodes.pdf(sx.tri_space)
1129 # -1=outside, 0=success, 1=failure, 2=mix
1130 args = {'event_id':2, 'indices':indices}
1131 sx.on_add_simplex(CandyNodes(nodes, args, pfv=pfv, pfw=pfw))
1132 else:
1133 fx = sx.f.sample_pdf(x.T, sx.tri_space)
1134
1135 return fx, fx*pfv, fx*pfw
1236
1237 x = quadpy.tn.transform(sx.tn_scheme.points, vertices_model.T)
1238 vol = quadpy.tn.get_vol(vertices_model)
1239
1240 nodes = sx.f.new_sample(x.T, sx.tri_space)
1241 fx = nodes.pdf(sx.tri_space)
1242
1243 p_mixed = vol * np.dot(fx, sx.tn_scheme.weights)
1136 1244
1137 1245
1138 simplex_measures = sx.tn_scheme.integrate(_get_pdf, vertices_model)
1139 #č místo fallback integrací
1140 mask = simplex_measures < 0
1141 if np.any(mask):
1142 print("Negative measures have occured in simplex %s" % indices)
1246 if p_mixed < 0:
1247 print("Negative measure has occured in simplex %s" % indices)
1143 1248 print("Rozbíla se nám integrace, integráční schema je na houby.") print("Rozbíla se nám integrace, integráční schema je na houby.")
1144 simplex_measures[mask] = 0
1249 p_mixed = 0
1145 1250
1146 p_mixed, pfv_simplex, pfw_simplex = simplex_measures
1147 1251
1252 pfv = p_mixed * fr
1253 pfw = p_mixed * wfr
1254 pf = vol * mfpdf
1148 1255
1149 1256 assert p_mixed >= 0 assert p_mixed >= 0
1150 assert pfv_simplex >= 0
1151 assert pfw_simplex >= 0
1257 assert pfv >= 0
1258 assert pfw >= 0
1259 assert pf >= 0
1152 1260
1153 1261 #č přídavám jednou, sčítám pořad dokola #č přídavám jednou, sčítám pořad dokola
1154 simplex_estimation = np.array((p_mixed, 0, pfv_simplex, pfw_simplex))
1155 sx.simplex_stats[tuple(indices)] = simplex_estimation
1262 simplex_estimation = np.array((p_mixed, pfv, pfw, pf))
1263 sx.mixed_simplices[simplex] = simplex_estimation
1264 #č odhady jsou ve slovníku, posíláme jen to, co tam není
1265 sx.on_mixed_added(simplex, indices, vertices_model, nodes, vol, fr, wfr, mfpdf)
1156 1266
1157 1267
1158 def _integrate_red_simplex(sx, indices):
1268 def _integrate_red_simplex(sx, simplex, indices):
1159 1269
1160 1270 vertices_model = sx.tri.points[indices] vertices_model = sx.tri.points[indices]
1161 1271
1162 # quadpy
1163 def _get_pdf(x):
1164
1165 # side effect
1166 if sx.on_add_simplex is not None:
1167 nodes = sx.f.new_sample(x.T, sx.tri_space)
1168 fx = nodes.pdf(sx.tri_space)
1169 # -1=outside, 0=success, 1=failure, 2=mix
1170 args = {'event_id':1, 'indices':indices}
1171 sx.on_add_simplex(CandyNodes(nodes, args))
1172 else:
1173 fx = sx.f.sample_pdf(x.T, sx.tri_space)
1174
1175 return fx
1272 x = quadpy.tn.transform(sx.tn_scheme.points, vertices_model.T)
1273 vol = quadpy.tn.get_vol(vertices_model)
1274
1275 nodes = sx.f.new_sample(x.T, sx.tri_space)
1276 fx = nodes.pdf(sx.tri_space)
1277
1278 p_failure = vol * np.dot(fx, sx.tn_scheme.weights)
1176 1279
1177 1280
1178 p_simplex = sx.tn_scheme.integrate(_get_pdf, vertices_model)
1179 #č místo fallback integrací
1180 if p_simplex < 0:
1181 print("Negative measures have occured in simplex %s" % indices)
1281 if p_failure < 0:
1282 print("Negative measure has occured in simplex %s" % indices)
1182 1283 print("Rozbíla se nám integrace, integráční schema je na houby.") print("Rozbíla se nám integrace, integráční schema je na houby.")
1183 p_simplex = 0
1284 p_failure = 0
1184 1285
1185 assert p_simplex >= 0
1186 #č přídavám jednou, sčítám pořad dokola
1187 simplex_estimation = np.array((0, p_simplex, p_simplex, p_simplex))
1188 sx.simplex_stats[tuple(indices)] = simplex_estimation
1286 assert p_failure >= 0
1287
1288 sx.failure_simplices[simplex] = p_failure
1289 #č odhady jsou ve slovníku, posíláme jen to, co tam není
1290 sx.on_failure_added(simplex, indices, vertices_model, nodes, vol)
1189 1291
1190 1292
1191 1293
1192 #č to je právý pf_estimation
1193 #č jiné třídy ve skutku vracej tri_estimation
1194 1294 def get_pf_estimation(sx): def get_pf_estimation(sx):
1195 # TRI-compatible estimation
1196 # -1=outside, 0=success, 1=failure, 2=mix
1197 # p_mixed, p_failure, p_fv, p_fw
1295 # p_mixed, pfv, pfw, pf
1198 1296 simplex_estimations = np.zeros(4, dtype=float) simplex_estimations = np.zeros(4, dtype=float)
1199 1297
1200 1298 #č nechce se mi počitat přes numpy: np.array(tuple(c.values())) #č nechce se mi počitat přes numpy: np.array(tuple(c.values()))
1201 1299 # let's iterate over all simplices we have # let's iterate over all simplices we have
1202 1300 # (calling code is responsible for the .simplex_stats validity) # (calling code is responsible for the .simplex_stats validity)
1203 for simplex_estimation in sx.simplex_stats.values():
1301 for simplex_estimation in sx.mixed_simplices.values():
1204 1302 simplex_estimations += simplex_estimation simplex_estimations += simplex_estimation
1205 1303
1206 p_mixed, p_failure, p_fv, p_fw = simplex_estimations
1207
1208 #ё так, для красоты
1209 global_stats = dict()
1210 # outside dodá Ghull
1211 global_stats['success_points'] = None #č další kód musí to přepsat
1212 global_stats['failure_points'] = None #č další kód musí to přepsat
1213 global_stats['success'] = None #č další kód musí to přepsat
1214 global_stats['failure'] = p_failure
1215 global_stats['mix'] = p_mixed
1216 global_stats['vertex_estimation'] = p_fv
1217 global_stats['weighted_vertex_estimation'] = p_fw
1218 global_stats['nsimplex'] = sx.tri.nsimplex
1219 global_stats['tn_scheme'] = sx.tn_scheme.name
1220 global_stats['tn_scheme_points'] = sx.tn_scheme.points.shape[1]
1221 global_stats['newly_invalidated'] = sx.newly_invalidated
1222 global_stats['newly_estimated'] = sx.newly_estimated
1223 global_stats['simplex_stats'] = len(sx.simplex_stats)
1224 global_stats['candidates_sets'] = None #č další kód musí to přepsat
1225 global_stats['ncoplanar'] = len(sx.tri.coplanar)
1226
1227 return global_stats
1304 p_failure = 0
1305 for failure_simplex in sx.failure_simplices.values():
1306 p_failure += failure_simplex
1307
1308
1309 p_mixed, pfv, pfw, pf = simplex_estimations
1310 nsim, nvar = sx.tri.points.shape
1311 return FastCubatureEstimation(
1312 nvar,
1313 nsim,
1314 p_failure,
1315 p_mixed,
1316 pfv + p_failure,
1317 pf + p_failure,
1318 pfw + p_failure,
1319 sx.tri.nsimplex,
1320 #sx.tn_scheme.name,
1321 sx.tn_scheme.points.shape[1],
1322 sx.newly_invalidated,
1323 sx.newly_estimated,
1324 len(sx.failure_simplices),
1325 len(sx.mixed_simplices),
1326 len(sx.tri.coplanar)
1327 )
1228 1328
1229 1329
1230 1330
Hints:
Before first commit, do not forget to setup your git environment:
git config --global user.name "your_name_here"
git config --global user.email "your@email_here"

Clone this repository using HTTP(S):
git clone https://rocketgit.com/user/iam-git/WellMet

Clone this repository using ssh (do not forget to upload a key first):
git clone ssh://rocketgit@ssh.rocketgit.com/user/iam-git/WellMet

Clone this repository using git:
git clone git://git.rocketgit.com/user/iam-git/WellMet

You are allowed to anonymously push to this repository.
This means that your pushed commits will automatically be transformed into a merge request:
... clone the repository ...
... make some changes and some commits ...
git push origin main