File wellmet/convex_hull.py changed (mode: 100644) (index c5f0c29..520661e) |
... |
... |
import mpmath |
6 |
6 |
from scipy import stats |
from scipy import stats |
7 |
7 |
from scipy import spatial # for QHull |
from scipy import spatial # for QHull |
8 |
8 |
|
|
|
9 |
|
from . import sball |
|
10 |
|
|
9 |
11 |
import quadpy |
import quadpy |
10 |
12 |
|
|
11 |
13 |
from collections import namedtuple |
from collections import namedtuple |
|
... |
... |
class CompleteHull: |
803 |
805 |
|
|
804 |
806 |
|
|
805 |
807 |
|
|
|
808 |
|
|
|
809 |
|
|
|
810 |
|
|
|
811 |
|
QHullEstimation = namedtuple('QHullEstimation', ( |
|
812 |
|
'nvar', |
|
813 |
|
'nsim', |
|
814 |
|
'nfacets', |
|
815 |
|
'r', 'R', |
|
816 |
|
'inner', |
|
817 |
|
'shell', |
|
818 |
|
'outer', |
|
819 |
|
'chi_outside', |
|
820 |
|
'orth_outside', |
|
821 |
|
#'tn_scheme', |
|
822 |
|
#'tn_scheme_points', |
|
823 |
|
'inside', |
|
824 |
|
'outside' |
|
825 |
|
)) |
|
826 |
|
|
|
827 |
|
|
806 |
828 |
class QHull: |
class QHull: |
807 |
829 |
# take some global function |
# take some global function |
808 |
830 |
_orth_helper = _orth_helper |
_orth_helper = _orth_helper |
|
... |
... |
class QHull: |
813 |
835 |
#get_2FORM_outside = get_2FORM_outside |
#get_2FORM_outside = get_2FORM_outside |
814 |
836 |
get_2FORM_equations = get_2FORM_equations |
get_2FORM_equations = get_2FORM_equations |
815 |
837 |
|
|
816 |
|
def __init__(self, sample, space='G', incremental=True, auto_update=True): |
|
|
838 |
|
def __init__(self, sample, space='G', incremental=True, auto_update=True, tn_scheme=None): |
817 |
839 |
self.sample = sample |
self.sample = sample |
818 |
840 |
self.incremental = incremental |
self.incremental = incremental |
819 |
841 |
self.space = space |
self.space = space |
820 |
842 |
self.auto_update = auto_update |
self.auto_update = auto_update |
|
843 |
|
self.tn_scheme = tn_scheme |
|
844 |
|
self.sball = sball.Sball(sample.nvar) |
|
845 |
|
self.shell = sball.Shell(sample.nvar) |
821 |
846 |
self.fallback_plan = quadpy.un.mysovskikh_1(sample.nvar).points |
self.fallback_plan = quadpy.un.mysovskikh_1(sample.nvar).points |
822 |
847 |
|
|
823 |
848 |
def regen(self): |
def regen(self): |
|
... |
... |
class QHull: |
965 |
990 |
else: |
else: |
966 |
991 |
return 0 |
return 0 |
967 |
992 |
|
|
|
993 |
|
def get_R(self): |
|
994 |
|
assert self.space == 'G' |
|
995 |
|
sum_squared = np.sum(np.square(self.sample.G), axis=1) |
|
996 |
|
#index = np.argmax(sum_squared) |
|
997 |
|
return np.sqrt(np.nanmax(sum_squared)) |
|
998 |
|
|
968 |
999 |
def shot(hull, ns, use_MC=False): |
def shot(hull, ns, use_MC=False): |
969 |
1000 |
try: |
try: |
970 |
1001 |
# take global function for shot() |
# take global function for shot() |
|
... |
... |
class QHull: |
997 |
1028 |
#č neřikejte mi, že musím pěčlivějc zpracovavat chyby! |
#č neřikejte mi, že musím pěčlivějc zpracovavat chyby! |
998 |
1029 |
return 1 |
return 1 |
999 |
1030 |
|
|
1000 |
|
def get_inside(hull): |
|
|
1031 |
|
def get_norm_inside(hull): |
|
1032 |
|
hull._update() |
1001 |
1033 |
if hull.space != 'G': |
if hull.space != 'G': |
1002 |
1034 |
raise |
raise |
1003 |
|
scheme = quadpy.tn.grundmann_moeller(hull.sample.nvar - 1, 5) |
|
1004 |
1035 |
|
|
1005 |
|
#np.sum(quadpy.tn.get_vol(box.G[box.convex_hull.simplices].transpose((1,0,2)))*np.sum(PQR * box.convex_hull.A, axis=1)) / 3 |
|
|
1036 |
|
vertices_model = hull.points[hull.simplices].transpose((1, 0, 2)) |
|
1037 |
|
result = hull.tn_scheme.integrate(hull._norm_inside_callback, vertices_model) |
1006 |
1038 |
|
|
|
1039 |
|
return np.sum(result, axis=1) |
1007 |
1040 |
|
|
1008 |
|
#np.sum(sheme.integrate(lambda x: np.full((26,1), 1), box.G[box.convex_hull.simplices].transpose((1, 0, 2)))) |
|
1009 |
|
#np.sum(sheme.integrate(hull._quadpy_callback, box.G[box.convex_hull.simplices].transpose((1, 0, 2))),axis=1) |
|
1010 |
1041 |
|
|
|
1042 |
|
def get_chi_outside(hull): |
|
1043 |
|
hull._update() |
|
1044 |
|
if hull.space != 'G': |
|
1045 |
|
raise |
|
1046 |
|
scheme = quadpy.tn.grundmann_moeller(hull.sample.nvar - 1, 5) |
1011 |
1047 |
|
|
1012 |
1048 |
vertices_model = hull.points[hull.simplices].transpose((1, 0, 2)) |
vertices_model = hull.points[hull.simplices].transpose((1, 0, 2)) |
|
1049 |
|
result = hull.tn_scheme.integrate(hull._chi_outside_callback, vertices_model) |
1013 |
1050 |
|
|
1014 |
|
return np.sum(scheme.integrate(hull._quadpy_callback, vertices_model), axis=1) |
|
|
1051 |
|
return np.sum(result) |
1015 |
1052 |
|
|
1016 |
1053 |
# quadpy |
# quadpy |
1017 |
|
def _quadpy_callback(hull, x): |
|
|
1054 |
|
def _norm_inside_callback(hull, x): |
1018 |
1055 |
# x.shape == (simplex_dim + 1, nsimplex, scheme_npoints) #(3, 26, 56) |
# x.shape == (simplex_dim + 1, nsimplex, scheme_npoints) #(3, 26, 56) |
1019 |
1056 |
|
|
1020 |
1057 |
cdfs = stats.norm.cdf(x) |
cdfs = stats.norm.cdf(x) |
|
... |
... |
class QHull: |
1028 |
1065 |
|
|
1029 |
1066 |
# n_values x nsimplex x scheme_points |
# n_values x nsimplex x scheme_points |
1030 |
1067 |
return cdfs |
return cdfs |
|
1068 |
|
|
|
1069 |
|
|
|
1070 |
|
# quadpy |
|
1071 |
|
def _chi_outside_callback(hull, x): |
|
1072 |
|
# x.shape == (simplex_dim + 1, nsimplex, scheme_npoints) #(3, 26, 56) |
|
1073 |
|
nvar, nsimplex, ninodes = x.shape |
1031 |
1074 |
|
|
|
1075 |
|
r = np.sqrt(np.sum(np.square(x), axis=0)) |
|
1076 |
|
cdfs = hull.sball.sf(r) |
|
1077 |
|
# chi_pdf * bůhvíco == phi(x) * phi(y) * phi(z) |
|
1078 |
|
chi_pdf = stats.chi.pdf(r, nvar) |
|
1079 |
|
norm_pdf = np.prod(stats.norm.pdf(x), axis=0) |
1032 |
1080 |
|
|
|
1081 |
|
# cdfs /= r |
|
1082 |
|
# cdfs /= nvar * 2 * np.pi |
1033 |
1083 |
|
|
|
1084 |
|
# n_values x nsimplex x scheme_points |
|
1085 |
|
return cdfs * norm_pdf / chi_pdf |
|
1086 |
|
|
|
1087 |
|
|
|
1088 |
|
|
|
1089 |
|
def get_convex_hull_estimation(hull): |
|
1090 |
|
|
|
1091 |
|
r = hull.get_r() |
|
1092 |
|
R = hull.get_R() |
|
1093 |
|
if r<0: |
|
1094 |
|
hull.shell.set_bounds(0, R) |
|
1095 |
|
else: |
|
1096 |
|
hull.shell.set_bounds(r, R) |
|
1097 |
|
|
|
1098 |
|
outer = hull.shell.pf |
|
1099 |
|
shell = hull.shell.p_shell |
|
1100 |
|
|
|
1101 |
|
chi_outside = hull.get_chi_outside() |
|
1102 |
|
|
|
1103 |
|
r_outside = outer + shell |
|
1104 |
|
if chi_outside > r_outside: |
|
1105 |
|
outside = r_outside |
|
1106 |
|
else: |
|
1107 |
|
outside = chi_outside |
|
1108 |
|
|
|
1109 |
|
return QHullEstimation( |
|
1110 |
|
hull.sample.nvar, |
|
1111 |
|
hull.npoints, |
|
1112 |
|
hull.nsimplex, |
|
1113 |
|
r, R, |
|
1114 |
|
hull.shell.ps, |
|
1115 |
|
shell, |
|
1116 |
|
outer, |
|
1117 |
|
chi_outside, |
|
1118 |
|
hull.get_orth_outside(), |
|
1119 |
|
#'tn_scheme', |
|
1120 |
|
#'tn_scheme_points', |
|
1121 |
|
1 - outside, |
|
1122 |
|
outside |
|
1123 |
|
) |
|
1124 |
|
|
|
1125 |
|
|
|
1126 |
|
|
|
1127 |
|
|
1034 |
1128 |
|
|
1035 |
1129 |
# Gaussian brick |
# Gaussian brick |
1036 |
1130 |
#č Gaussovská cihla |
#č Gaussovská cihla |
|
... |
... |
class QHullCubature(QHull): |
1314 |
1408 |
else: |
else: |
1315 |
1409 |
self.regen() |
self.regen() |
1316 |
1410 |
|
|
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 |
1411 |
|
|
1323 |
1412 |
def get_inside(hull): |
def get_inside(hull): |
1324 |
1413 |
hull._update() |
hull._update() |