File wellmet/convex_hull.py changed (mode: 100644) (index 340a019..c5f0c29) |
... |
... |
from scipy import spatial # for QHull |
8 |
8 |
|
|
9 |
9 |
import quadpy |
import quadpy |
10 |
10 |
|
|
|
11 |
|
from collections import namedtuple |
|
12 |
|
|
11 |
13 |
mpmath.mp.dps = 325 # to cover anything that double-presigion float can handle |
mpmath.mp.dps = 325 # to cover anything that double-presigion float can handle |
12 |
14 |
|
|
13 |
15 |
# maximum radius, where norm.pdf() wasn't zero |
# maximum radius, where norm.pdf() wasn't zero |
|
... |
... |
class Grick: |
1262 |
1264 |
return np.min((-hull.bm, hull.bp)) |
return np.min((-hull.bm, hull.bp)) |
1263 |
1265 |
|
|
1264 |
1266 |
|
|
|
1267 |
|
|
|
1268 |
|
|
|
1269 |
|
|
|
1270 |
|
|
|
1271 |
|
QHullCubatureEstimation = namedtuple('QHullCubatureEstimation', ( |
|
1272 |
|
'nvar', |
|
1273 |
|
'nsim', |
|
1274 |
|
'r', 'R', |
|
1275 |
|
'nfacets', |
|
1276 |
|
#'tn_scheme', |
|
1277 |
|
'tn_scheme_points', |
|
1278 |
|
'inside', |
|
1279 |
|
'outside' |
|
1280 |
|
)) |
|
1281 |
|
|
|
1282 |
|
|
|
1283 |
|
class QHullCubature(QHull): |
|
1284 |
|
|
|
1285 |
|
def __init__(hull, sample, tn_scheme, space='G', incremental=True, auto_update=True): |
|
1286 |
|
|
|
1287 |
|
# for now |
|
1288 |
|
#č jednotka dočasnosti - jeden furt |
|
1289 |
|
assert space == 'G' |
|
1290 |
|
|
|
1291 |
|
hull.tn_scheme = tn_scheme |
|
1292 |
|
super().__init__(sample, space, incremental, auto_update) |
|
1293 |
|
|
|
1294 |
|
|
|
1295 |
|
def regen(self): |
|
1296 |
|
points = getattr(self.sample, self.space) |
|
1297 |
|
self.convex_hull = spatial.ConvexHull(points, incremental=self.incremental) |
|
1298 |
|
self.facets = dict() |
|
1299 |
|
|
|
1300 |
|
self.newly_estimated = 0 |
|
1301 |
|
self.newly_invalidated = 0 |
|
1302 |
|
|
|
1303 |
|
simplices = self.convex_hull.simplices |
|
1304 |
|
x_normals = self.convex_hull.equations[:,0] |
|
1305 |
|
|
|
1306 |
|
for indices, normal_x in zip(simplices, x_normals): |
|
1307 |
|
facet = tuple(indices) |
|
1308 |
|
self._integrate_facet(facet, indices, normal_x) |
|
1309 |
|
|
|
1310 |
|
def update(self): |
|
1311 |
|
if self.enough_points and (self.convex_hull.npoints < self.sample.nsim): |
|
1312 |
|
if self.incremental: |
|
1313 |
|
self.__update() |
|
1314 |
|
else: |
|
1315 |
|
self.regen() |
|
1316 |
|
|
|
1317 |
|
def get_R(self): |
|
1318 |
|
sum_squared = np.sum(np.square(self.sample.G), axis=1) |
|
1319 |
|
#index = np.argmax(sum_squared) |
|
1320 |
|
return np.sqrt(np.nanmax(sum_squared)) |
|
1321 |
|
|
|
1322 |
|
|
|
1323 |
|
def get_inside(hull): |
|
1324 |
|
hull._update() |
|
1325 |
|
if hull.enough_points: |
|
1326 |
|
p_inside = mpmath.mpf(0) |
|
1327 |
|
for p_facet in hull.facets.values(): |
|
1328 |
|
p_inside += p_facet |
|
1329 |
|
return p_inside |
|
1330 |
|
else: |
|
1331 |
|
return 0 |
|
1332 |
|
|
|
1333 |
|
|
|
1334 |
|
|
|
1335 |
|
def __update(hull): |
|
1336 |
|
|
|
1337 |
|
hull._hull_update() #č jako vždy, chyby nechytáme |
|
1338 |
|
|
|
1339 |
|
#č podstata problému byla a je v tom, |
|
1340 |
|
#č že numpy neumí pracovat s věktory, případně s submaticemi |
|
1341 |
|
#č jako s celky. Má routiny pouze pro 1D množiny. |
|
1342 |
|
#č Navíc, numpy vektory nejsou hashable, nejde je přímo hodit |
|
1343 |
|
#č ani do setů, ani do slovníků |
|
1344 |
|
#č Nezbyvá nám nic jineho, než prévest ndarray na množinu n-tic. |
|
1345 |
|
|
|
1346 |
|
#č na rozdil od triangulace, zde potřebujeme i normály, aspoň první složku |
|
1347 |
|
#č jenomže výtahnout ji z numpy equations později na zakladě indexů by bylo velmi dráho. |
|
1348 |
|
new_facets = dict(zip((tuple(simplex) for simplex in hull.simplices), hull.equations[:,0])) |
|
1349 |
|
|
|
1350 |
|
# possible naive implementation |
|
1351 |
|
#to_invalidate = sx.simplices_set - new_simplices_set |
|
1352 |
|
#to_estimate = new_simplices_set - sx.simplices_set |
|
1353 |
|
|
|
1354 |
|
facets = hull.facets # dict |
|
1355 |
|
#č vymejšlím, jak by se mohlo ušetřit, |
|
1356 |
|
#č aby se neprobíhala smyčka přes obrovský set na dvakrat |
|
1357 |
|
# Update the set, keeping only elements found in either set, but not in both |
|
1358 |
|
either_facets = new_facets.keys() ^ facets.keys() # set |
|
1359 |
|
to_estimate = either_facets - facets.keys() # set |
|
1360 |
|
to_invalidate = either_facets - to_estimate # set |
|
1361 |
|
|
|
1362 |
|
hull.newly_estimated = len(to_estimate) |
|
1363 |
|
hull.newly_invalidated = len(to_invalidate) |
|
1364 |
|
|
|
1365 |
|
#č invalidace |
|
1366 |
|
for key in to_invalidate: |
|
1367 |
|
del facets[key] |
|
1368 |
|
|
|
1369 |
|
|
|
1370 |
|
for facet in to_estimate: |
|
1371 |
|
indices = np.array(facet) |
|
1372 |
|
normal_x = new_facets[facet] |
|
1373 |
|
hull._integrate_facet(facet, indices, normal_x) |
|
1374 |
|
|
|
1375 |
|
|
|
1376 |
|
|
|
1377 |
|
def _hull_update(self): |
|
1378 |
|
points = getattr(self.sample, self.space) |
|
1379 |
|
self.convex_hull.add_points(points[self.convex_hull.npoints:]) |
|
1380 |
|
|
|
1381 |
|
|
|
1382 |
|
|
|
1383 |
|
def get_convex_hull_estimation(hull): |
|
1384 |
|
inside = hull.get_inside() |
|
1385 |
|
return QHullCubatureEstimation( |
|
1386 |
|
hull.sample.nvar, |
|
1387 |
|
hull.npoints, |
|
1388 |
|
hull.get_r(), hull.get_R(), |
|
1389 |
|
hull.nsimplex, |
|
1390 |
|
#'tn_scheme', |
|
1391 |
|
hull.tn_scheme.points.shape[1], |
|
1392 |
|
float(inside), |
|
1393 |
|
float(1 - inside) |
|
1394 |
|
) |
1265 |
1395 |
|
|
|
1396 |
|
|
|
1397 |
|
def _integrate_facet(hull, facet, indices, normal_x): |
|
1398 |
|
|
|
1399 |
|
|
|
1400 |
|
vertices_model = hull.points[indices] |
1266 |
1401 |
|
|
|
1402 |
|
# nvar x n_integration_nodes |
|
1403 |
|
x = quadpy.tn.transform(hull.tn_scheme.points, vertices_model.T) |
|
1404 |
|
vol = quadpy.tn.get_vol(vertices_model) |
|
1405 |
|
|
|
1406 |
|
if not np.isfinite(vol): |
|
1407 |
|
print("Incorrect volume has occured in simplex %s" % indices) |
|
1408 |
|
hull.facets[facet] = 0 |
|
1409 |
|
return |
|
1410 |
|
|
|
1411 |
|
nvar = hull.sample.nvar |
|
1412 |
|
store = mpmath.mpf(0) |
|
1413 |
|
|
|
1414 |
|
for node, weight in zip(x.T, hull.tn_scheme.weights): |
|
1415 |
|
node_iterator = iter(node) |
|
1416 |
|
node_store = mpmath.ncdf(next(node_iterator)) |
|
1417 |
|
for __ in range(nvar - 1): |
|
1418 |
|
node_store *= mpmath.npdf(next(node_iterator)) |
|
1419 |
|
|
|
1420 |
|
store += node_store * weight |
|
1421 |
|
|
|
1422 |
|
|
|
1423 |
|
if store < 0: |
|
1424 |
|
print("Negative measure has occured in facet %s" % indices) |
|
1425 |
|
hull.facets[facet] = 0 |
|
1426 |
|
return |
|
1427 |
|
|
|
1428 |
|
hull.facets[facet] = store * normal_x * vol |
|
1429 |
|
|
|
1430 |
|
|
|
1431 |
|
|
|
1432 |
|
|
|
1433 |
|
|
|
1434 |
|
|
|
1435 |
|
|
|
1436 |
|
|
|
1437 |
|
|
|
1438 |
|
|
|
1439 |
|
|
1267 |
1440 |
|
|