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)
convex_hull.QHullCubature: implement memory-friendly outside (chi) integration ad8210a04b1e0903de7435cad16b1304707d0e6e I am 2023-02-09 22:22:05
qt_gui.qt_plot: require box recommendation to show next point 6f726047f7f08e884359020eaa1eac6f6cc125d2 I am 2023-02-09 11:51:44
dicebox.circumtri.CirQTri.get_circum_node: limit circumcenter by explore radius, not by just R 136ec73bb06da16c1f2bce64b3c349be4c8ba975 I am 2023-02-09 03:09:51
dicebox.circumtri: implement new brand CirQTri box 5879b8ad6317c65aa9481b59f76b6159f19e247a I am 2023-02-09 01:29:10
simplex.FullCubatureIntegration: store simplex probabilities in sorted dicts c0da90731ff3ede47d9b4eec0ad9b28a29027167 I am 2023-02-09 01:23:14
dicebox.circumtri: exploratory: even better idea 811ab11cd7172ff4a3807992f4928be2e8068ec0 I am 2023-02-08 15:31:23
dicebox.circumtri: exploratory, new idea 526d3f6185887ff48a341b0705d74dde4d15ca87 I am 2023-02-08 03:03:41
dicebox.circumtri: exploratory 806063d2e223c812280dc4845153450dd47faed3 I am 2023-02-06 17:15:15
dicebox.circumtri: exploratory efed2589f642d502f30e80f0e9b45dfeecd1c7c7 I am 2023-02-06 13:40:24
dicebox.circumtri: exploratory - again 34d3f4e47420e1c1e26b09570fb44d3270194751 I am 2023-02-06 12:50:45
qt_gui.qt_dicebox: change default q of circumtri classes 9fd5855e5d7cacf80d27fb383dd18a92d60e138b I am 2023-02-06 12:30:27
dicebox.circumtri: tune up exploratory rule one more time bfaa8d65bd13956252a6a25382c621aca7a33e3f I am 2023-02-06 12:05:00
convex_hull.QHull.is_outside: mark everything as outside until convex hull is created beecd924a63603868807a8a603cbeb04857e5cde I am 2023-02-06 12:03:35
dicebox.circumtri: get_exploratory_radius: turn it up some more ecd638ca90bac88df515cf803a103e350b08041e I am 2023-02-03 10:58:37
mplot: ongoing ad-hoc changes fd71b93abee959231ffd76324cf2fbc0b3bbbdfe I am 2023-02-02 17:59:27
dicebox.circumtri: change default q setting 2a9bdc6604113462a935d409bca429b49ac52cbc I am 2023-02-02 11:55:41
dicebox.circumtri: minor change in explore radia 9aa06dc11be7c935dd47715c0aa53b401be951d4 I am 2023-02-01 23:37:13
dicebox.circumtri: exploratory: another iteration 78c849a2c1b07a55f87d32e3cd6570db222a50b6 I am 2023-02-01 01:02:00
dicebox: new experimental explore() 1d570ef1dcc66b9b1015bd6c1f6877e1405e676c I am 2023-02-01 00:34:57
mplot: improve GRaph; draw vertex estimator only in convergence diagram 5c22e808e4c6cabd8e0d077d5f295eea93957582 I am 2023-01-31 18:02:00
Commit ad8210a04b1e0903de7435cad16b1304707d0e6e - convex_hull.QHullCubature: implement memory-friendly outside (chi) integration
Author: I am
Author date (UTC): 2023-02-09 22:22
Committer name: I am
Committer date (UTC): 2023-02-09 22:22
Parent(s): 6f726047f7f08e884359020eaa1eac6f6cc125d2
Signer:
Signing key:
Signing status: N
Tree: cebc01e7ac1e77aef884a2d1921ebd500b9c745b
File Lines added Lines deleted
wellmet/convex_hull.py 81 162
wellmet/dicebox/circumtri.py 2 2
File wellmet/convex_hull.py changed (mode: 100644) (index a48c4d2..e8774b9)
... ... class QHull:
1221 1221 ) )
1222 1222
1223 1223
1224 class QHullCubature(QHull):
1225
1226 def get_chi_outside(hull):
1227 hull._update()
1228
1229 if hull.space != 'G':
1230 raise
1231
1232
1233 integrals = 0
1234 # if first time runned
1235 if 'facets' not in dir(hull):
1236 hull.facets = facets = dict()
1237 ncubature = len(hull.tn_scheme.weights)
1238 hull._x = np.empty((hull.sample.nvar, ncubature))
1239 hull._r = np.empty(ncubature)
1240 hull._result = np.empty(ncubature)
1241 for indices, normal in zip(hull.simplices, hull.A):
1242 facet_integral = hull.get_facet_outside(indices, normal)
1243 facets[indices.tobytes()] = facet_integral
1244 integrals += facet_integral
1245
1246 return integrals / hull.nsphere_surface_area
1247
1248 facets = hull.facets
1249 new_facet_dict = dict()
1250 for indices, normal in zip(hull.simplices, hull.A):
1251 key = indices.tobytes()
1252 if key in facets:
1253 facet_integral = facets[key]
1254 else:
1255 facet_integral = hull.get_facet_outside(indices, normal)
1256 new_facet_dict[key] = facet_integral
1257 integrals += facet_integral
1258
1259 hull.facets = new_facet_dict
1260 return integrals / hull.nsphere_surface_area
1261
1262
1263
1264
1265 def get_facet_outside(hull, indices, normal):
1266 vertices_model = hull.points[indices]
1267
1268 vol = quadpy.tn.get_vol(vertices_model)
1269 if not np.isfinite(vol):
1270 print("QHull: Incorrect area has occured in facet %s" % indices)
1271 return 0
1272
1273 nvar = len(normal)
1274 x = hull._x
1275 r_nvar = hull._r
1276 result = hull._result
1277
1278 # nvar x n_integration_nodes
1279 # Transform the points `xi` from the reference simplex onto `simplex`.
1280 # same as quadpy.tn.transform(hull.tn_scheme.points, vertices_model.T)
1281 #np.dot(simplex, points)
1282 np.matmul(vertices_model.T, hull.tn_scheme.points, out=x)
1283
1284 np.einsum('ij,ij->j', x, x, out=result) # r**2
1285 np.power(result, nvar/2, out=r_nvar) # r**(nvar-1) (chi) * r (to normalize)
1286
1287 np.divide(result, 2, out=result) # r**2 / 2
1288 special.gammaincc(nvar / 2, result, out=result) # F
1289
1290 #č normalizujeme normály. Jenže normály samotné dodáme později
1291 result /= r_nvar
1292
1293 #č normála
1294 np.dot(normal, x, out=r_nvar)
1295 result *= r_nvar
1296
1297 # facet_outside = integral / unit_ball_surface_area
1298 integral = vol * np.inner(result, hull.tn_scheme.weights)
1299 if integral > 0:
1300 return integral
1301 else:
1302 print("QHull: Negative value has appeared in facet %s during integration" % indices)
1303 return 0
1224 1304
1225 1305
1226 1306
 
... ... class Grick:
1460 1540
1461 1541
1462 1542
1463 QHullCubatureEstimation = namedtuple('QHullCubatureEstimation', (
1464 'nvar',
1465 'nsim',
1466 'r', 'R',
1467 'nfacets',
1468 #'tn_scheme',
1469 'tn_scheme_points',
1470 'inside',
1471 'outside'
1472 ))
1473
1474
1475 class QHullCubature(QHull):
1476
1477 def __init__(hull, sample, tn_scheme, space='G', incremental=True, auto_update=True):
1478
1479 # for now
1480 #č jednotka dočasnosti - jeden furt
1481 assert space == 'G'
1482
1483 hull.tn_scheme = tn_scheme
1484 super().__init__(sample, space, incremental, auto_update)
1485
1486
1487 def regen(self):
1488 points = getattr(self.sample, self.space)
1489 self.convex_hull = spatial.ConvexHull(points, incremental=self.incremental)
1490 self.facets = dict()
1491
1492 self.newly_estimated = 0
1493 self.newly_invalidated = 0
1494
1495 simplices = self.convex_hull.simplices
1496 x_normals = self.convex_hull.equations[:,0]
1497
1498 for indices, normal_x in zip(simplices, x_normals):
1499 facet = tuple(indices)
1500 self._integrate_facet(facet, indices, normal_x)
1501
1502 def update(self):
1503 if self.enough_points and (self.convex_hull.npoints < self.sample.nsim):
1504 if self.incremental:
1505 self.__update()
1506 else:
1507 self.regen()
1508
1509
1510 def get_inside(hull):
1511 hull._update()
1512 if hull.enough_points:
1513 p_inside = mpmath.mpf(0)
1514 for p_facet in hull.facets.values():
1515 p_inside += p_facet
1516 return p_inside
1517 else:
1518 return 0
1519
1520
1521
1522 def __update(hull):
1523
1524 hull._hull_update() #č jako vždy, chyby nechytáme
1525
1526 #č podstata problému byla a je v tom,
1527 #č že numpy neumí pracovat s věktory, případně s submaticemi
1528 #č jako s celky. Má routiny pouze pro 1D množiny.
1529 #č Navíc, numpy vektory nejsou hashable, nejde je přímo hodit
1530 #č ani do setů, ani do slovníků
1531 #č Nezbyvá nám nic jineho, než prévest ndarray na množinu n-tic.
1532
1533 #č na rozdil od triangulace, zde potřebujeme i normály, aspoň první složku
1534 #č jenomže výtahnout ji z numpy equations později na zakladě indexů by bylo velmi dráho.
1535 new_facets = dict(zip((tuple(simplex) for simplex in hull.simplices), hull.equations[:,0]))
1536
1537 # possible naive implementation
1538 #to_invalidate = sx.simplices_set - new_simplices_set
1539 #to_estimate = new_simplices_set - sx.simplices_set
1540
1541 facets = hull.facets # dict
1542 #č vymejšlím, jak by se mohlo ušetřit,
1543 #č aby se neprobíhala smyčka přes obrovský set na dvakrat
1544 # Update the set, keeping only elements found in either set, but not in both
1545 either_facets = new_facets.keys() ^ facets.keys() # set
1546 to_estimate = either_facets - facets.keys() # set
1547 to_invalidate = either_facets - to_estimate # set
1548
1549 hull.newly_estimated = len(to_estimate)
1550 hull.newly_invalidated = len(to_invalidate)
1551
1552 #č invalidace
1553 for key in to_invalidate:
1554 del facets[key]
1555
1556
1557 for facet in to_estimate:
1558 indices = np.array(facet)
1559 normal_x = new_facets[facet]
1560 hull._integrate_facet(facet, indices, normal_x)
1561
1562
1563
1564 def _hull_update(self):
1565 points = getattr(self.sample, self.space)
1566 self.convex_hull.add_points(points[self.convex_hull.npoints:])
1567
1568
1569
1570 def get_convex_hull_estimation(hull):
1571 inside = hull.get_inside()
1572 return QHullCubatureEstimation(
1573 hull.sample.nvar,
1574 hull.npoints,
1575 hull.get_r(), hull.get_R(),
1576 hull.nsimplex,
1577 #'tn_scheme',
1578 hull.tn_scheme.points.shape[1],
1579 float(inside),
1580 float(1 - inside)
1581 )
1582
1583
1584 def _integrate_facet(hull, facet, indices, normal_x):
1585
1586
1587 vertices_model = hull.points[indices]
1588
1589 # nvar x n_integration_nodes
1590 x = quadpy.tn.transform(hull.tn_scheme.points, vertices_model.T)
1591 vol = quadpy.tn.get_vol(vertices_model)
1592
1593 if not np.isfinite(vol):
1594 print("Incorrect volume has occured in simplex %s" % indices)
1595 hull.facets[facet] = 0
1596 return
1597
1598 nvar = hull.sample.nvar
1599 store = mpmath.mpf(0)
1600
1601 for node, weight in zip(x.T, hull.tn_scheme.weights):
1602 node_iterator = iter(node)
1603 node_store = mpmath.ncdf(next(node_iterator))
1604 for __ in range(nvar - 1):
1605 node_store *= mpmath.npdf(next(node_iterator))
1606
1607 store += node_store * weight
1608
1609
1610 if store < 0:
1611 print("Negative measure has occured in facet %s" % indices)
1612 hull.facets[facet] = 0
1613 return
1614
1615 hull.facets[facet] = store * normal_x * vol
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625 1543
1544
1626 1545
1627 1546
File wellmet/dicebox/circumtri.py changed (mode: 100644) (index d633298..f357382)
... ... class CirQTri(_CircumTri):
378 378
379 379 def _regen_outside(bx): def _regen_outside(bx):
380 380 facet_scheme = quadpy.tn.grundmann_moeller(bx.nvar - 1, bx.convex_hull_degree) facet_scheme = quadpy.tn.grundmann_moeller(bx.nvar - 1, bx.convex_hull_degree)
381 bx.convex_hull = khull.QHull(bx.f_model, space='G', incremental=True,
381 bx.convex_hull = khull.QHullCubature(bx.f_model, space='G', incremental=True,
382 382 auto_update=False, tn_scheme=facet_scheme) auto_update=False, tn_scheme=facet_scheme)
383 383
384 384 #č konečně mám pořádnou stejtful třídu #č konečně mám pořádnou stejtful třídu
 
... ... class QTri(_CircumTri, _PsiTri):
940 940
941 941 def _regen_outside(bx): def _regen_outside(bx):
942 942 facet_scheme = quadpy.tn.grundmann_moeller(bx.nvar - 1, bx.convex_hull_degree) facet_scheme = quadpy.tn.grundmann_moeller(bx.nvar - 1, bx.convex_hull_degree)
943 bx.convex_hull = khull.QHull(bx.f_model, space='G', incremental=True,
943 bx.convex_hull = khull.QHullCubature(bx.f_model, space='G', incremental=True,
944 944 auto_update=False, tn_scheme=facet_scheme) auto_update=False, tn_scheme=facet_scheme)
945 945
946 946 #č konečně mám pořádnou stejtful třídu #č konečně mám pořádnou stejtful třídu
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