File wellmet/convex_hull.py changed (mode: 100644) (index 066a2d1..fae1de0) |
... |
... |
def _get_nD_orth_basis(hull): |
173 |
173 |
#č bůhví, jestli obálka negeneruje A a b-čko dynamicky |
#č bůhví, jestli obálka negeneruje A a b-čko dynamicky |
174 |
174 |
b = hull.b |
b = hull.b |
175 |
175 |
A = hull.A |
A = hull.A |
176 |
|
ndim = hull.sample.nvar |
|
|
176 |
|
|
|
177 |
|
return _get_nD_orth_basis_from_normals(A, b) |
|
178 |
|
|
|
179 |
|
|
|
180 |
|
def _get_nD_orth_basis_from_normals(A, b): |
|
181 |
|
_nplates, ndim = A.shape |
177 |
182 |
|
|
178 |
183 |
#č pro 2D "nastav první a pak po veškerem zbytku |
#č pro 2D "nastav první a pak po veškerem zbytku |
179 |
184 |
#č ještě ten předposlední vektor" špatně funguje |
#č ještě ten předposlední vektor" špatně funguje |
|
... |
... |
def get_2FORM_equations(hull): |
313 |
318 |
b = np.array([[x_backward],[-x_forward]]) |
b = np.array([[x_backward],[-x_forward]]) |
314 |
319 |
return np.hstack((A,b)) |
return np.hstack((A,b)) |
315 |
320 |
|
|
|
321 |
|
def get_bp_bm(A, points): |
|
322 |
|
B = A @ points.T |
|
323 |
|
bp = np.nanmax(B, axis=1) |
|
324 |
|
bm = np.nanmin(B, axis=1) |
|
325 |
|
return bp, bm |
|
326 |
|
|
316 |
327 |
|
|
317 |
328 |
#č jistě musíme mít nějaký zbytečný kus kódu |
#č jistě musíme mít nějaký zbytečný kus kódu |
318 |
329 |
#č třida jen pro formu, jen tak na hračku |
#č třida jen pro formu, jen tak na hračku |
|
... |
... |
class QHull: |
918 |
929 |
return 1 |
return 1 |
919 |
930 |
|
|
920 |
931 |
|
|
|
932 |
|
# Gaussian brick |
|
933 |
|
#č Gaussovská cihla |
|
934 |
|
# |
|
935 |
|
# it was a mistake to design previous classes without explicit update() |
|
936 |
|
class Grick: |
|
937 |
|
# take some global functions |
|
938 |
|
fire = fire |
|
939 |
|
shot = shot |
|
940 |
|
|
|
941 |
|
|
|
942 |
|
def get_FORM_outside(self): |
|
943 |
|
return stats.norm.sf(self.get_r()) |
|
944 |
|
|
|
945 |
|
|
|
946 |
|
def __init__(hull, sample, direct_plan, nrandom=50, auto_update=True): |
|
947 |
|
hull.sample = sample |
|
948 |
|
hull.direct_plan = direct_plan |
|
949 |
|
hull.space = 'G' |
|
950 |
|
hull.ndim = sample.nvar |
|
951 |
|
hull.nrandom = nrandom #č počet náhodných normalových věktorů |
|
952 |
|
|
|
953 |
|
hull.auto_update = auto_update |
|
954 |
|
|
|
955 |
|
hull.regen() |
|
956 |
|
|
|
957 |
|
|
|
958 |
|
def regen(hull): |
|
959 |
|
hull._npoints = 0 |
|
960 |
|
hull.mins = np.full(hull.ndim, np.inf) |
|
961 |
|
hull.maxs = np.full(hull.ndim, -np.inf) |
|
962 |
|
|
|
963 |
|
hull.bp_dp = np.full(len(hull.direct_plan), -np.inf) |
|
964 |
|
hull.bm_dp = np.full(len(hull.direct_plan), np.inf) |
|
965 |
|
|
|
966 |
|
hull.bp = np.full(hull.ndim, -np.inf) |
|
967 |
|
hull.bm = np.full(hull.ndim, np.inf) |
|
968 |
|
|
|
969 |
|
hull.orth_basis = np.eye(hull.ndim) |
|
970 |
|
|
|
971 |
|
hull.update() |
|
972 |
|
|
|
973 |
|
|
|
974 |
|
def _update(hull): |
|
975 |
|
if hull.auto_update: |
|
976 |
|
hull.update() |
|
977 |
|
|
|
978 |
|
# it was a mistake to design previous classes without explicit update() |
|
979 |
|
def update(hull): |
|
980 |
|
if hull._npoints < hull.npoints: |
|
981 |
|
hull._points = hull.sample.G |
|
982 |
|
|
|
983 |
|
new_points = hull._points[hull._npoints:] |
|
984 |
|
|
|
985 |
|
new_mins = np.nanmin(new_points, axis=0) |
|
986 |
|
new_maxs = np.nanmax(new_points, axis=0) |
|
987 |
|
|
|
988 |
|
np.nanmin((new_mins, hull.mins), axis=0, out=hull.mins) |
|
989 |
|
np.nanmax((new_maxs, hull.maxs), axis=0, out=hull.maxs) |
|
990 |
|
|
|
991 |
|
new_bp, new_bm = get_bp_bm(hull.direct_plan, new_points) |
|
992 |
|
np.nanmax((new_bp, hull.bp_dp), axis=0, out=hull.bp_dp) |
|
993 |
|
np.nanmin((new_bm, hull.bm_dp), axis=0, out=hull.bm_dp) |
|
994 |
|
|
|
995 |
|
|
|
996 |
|
new_bp, new_bm = get_bp_bm(hull.orth_basis, new_points) |
|
997 |
|
np.nanmax((new_bp, hull.bp), axis=0, out=hull.bp) |
|
998 |
|
np.nanmin((new_bm, hull.bm), axis=0, out=hull.bm) |
|
999 |
|
|
|
1000 |
|
hull._gen_orth_basis() |
|
1001 |
|
|
|
1002 |
|
hull._npoints = hull.npoints |
|
1003 |
|
|
|
1004 |
|
|
|
1005 |
|
def gen_random_Ab(hull): |
|
1006 |
|
ns = hull.ndim * 2 + len(hull.direct_plan) + hull.nrandom |
|
1007 |
|
#č necháme vygenerovat maticu náhodných čisel |
|
1008 |
|
#č ale použijeme jen část z ní |
|
1009 |
|
#č moje zkušenost je, že np.ones a np.random.randn běží stejně rychle |
|
1010 |
|
A = np.random.randn(ns, hull.ndim) #random directions |
|
1011 |
|
A[:hull.ndim] = np.eye(hull.ndim) |
|
1012 |
|
A[hull.ndim:hull.ndim*2] = hull.orth_basis |
|
1013 |
|
A[hull.ndim*2:ns-hull.nrandom] = hull.direct_plan |
|
1014 |
|
|
|
1015 |
|
# rand_dir: prepare ns random directions on a unit d-sphere |
|
1016 |
|
rand_dir = A[ns-hull.nrandom:] |
|
1017 |
|
|
|
1018 |
|
lengths = np.sum(np.square(rand_dir), axis=1) |
|
1019 |
|
lengths = np.sqrt(lengths, out=lengths) #lengths of each radius-vector |
|
1020 |
|
|
|
1021 |
|
# scale all radii-vectors to unit length |
|
1022 |
|
# use [:,None] to get an transposed 2D array |
|
1023 |
|
np.divide(rand_dir, lengths[:,None], out=rand_dir) |
|
1024 |
|
|
|
1025 |
|
b = np.empty(ns) |
|
1026 |
|
#č Nám jde o basis. Dáme ty nejlepší b-čka. |
|
1027 |
|
#č znaménka vektorů, tj. směry dopředu nebo dozadu v této fázi nevadí |
|
1028 |
|
np.nanmax((-hull.maxs, hull.mins), axis=0, out=b[:hull.ndim]) |
|
1029 |
|
np.nanmax((-hull.bp, hull.bm), axis=0, out=b[hull.ndim:hull.ndim*2]) |
|
1030 |
|
np.nanmax((-hull.bp_dp, hull.bm_dp), axis=0, out=b[hull.ndim*2:ns-hull.nrandom]) |
|
1031 |
|
|
|
1032 |
|
rbp, rbm = get_bp_bm(rand_dir, hull._points) |
|
1033 |
|
np.nanmax((-rbp, rbm), axis=0, out=b[ns-hull.nrandom:]) |
|
1034 |
|
|
|
1035 |
|
return A, b |
|
1036 |
|
|
|
1037 |
|
|
|
1038 |
|
def _gen_orth_basis(hull): |
|
1039 |
|
A, b = hull.gen_random_Ab() |
|
1040 |
|
|
|
1041 |
|
#č je to špatné z hlediska testování, |
|
1042 |
|
#č ale nevím jak z toho. |
|
1043 |
|
#č pro 2D musím udělat vyjimku z "předposlední" |
|
1044 |
|
if hull.ndim > 2: |
|
1045 |
|
hull.orth_basis = _get_nD_orth_basis_from_normals(A, b) |
|
1046 |
|
else: |
|
1047 |
|
to_fire = np.nanargmax(b) |
|
1048 |
|
a = A[to_fire] |
|
1049 |
|
x, y = a |
|
1050 |
|
hull.orth_basis = np.array([[x, y],[-y, x]]) |
|
1051 |
|
|
|
1052 |
|
|
|
1053 |
|
#č musí tam bejt G coordinates |
|
1054 |
|
hull.bp, hull.bm = get_bp_bm(hull.orth_basis, hull._points) |
|
1055 |
|
|
|
1056 |
|
|
|
1057 |
|
|
|
1058 |
|
|
|
1059 |
|
|
|
1060 |
|
|
|
1061 |
|
|
|
1062 |
|
|
|
1063 |
|
def get_orth_outside(hull): |
|
1064 |
|
hull._update() |
|
1065 |
|
return calculate_brick_complement_probability(hull.bm, hull.bp) |
|
1066 |
|
|
|
1067 |
|
|
|
1068 |
|
def get_orth_equations(hull): |
|
1069 |
|
hull._update() |
|
1070 |
|
return hull.equations |
|
1071 |
|
|
|
1072 |
|
|
|
1073 |
|
def get_2FORM_outside(hull): |
|
1074 |
|
hull._update() |
|
1075 |
|
return stats.norm.cdf(hull.bm[0]) + stats.norm.sf(hull.bp[0]) |
|
1076 |
|
|
|
1077 |
|
def get_2FORM_equations(hull): |
|
1078 |
|
hull._update() |
|
1079 |
|
A = np.vstack((-hull.orth_basis[0], hull.orth_basis[0])) |
|
1080 |
|
b = np.array([[hull.bm[0]],[-hull.bp[0]]]) |
|
1081 |
|
return np.hstack((A,b)) |
|
1082 |
|
|
|
1083 |
|
|
|
1084 |
|
@property |
|
1085 |
|
def points(hull): |
|
1086 |
|
if hull.auto_update: |
|
1087 |
|
return getattr(hull.sample, hull.space) |
|
1088 |
|
else: |
|
1089 |
|
return hull._points |
|
1090 |
|
|
|
1091 |
|
|
|
1092 |
|
@property |
|
1093 |
|
def npoints(hull): |
|
1094 |
|
if hull.auto_update: |
|
1095 |
|
return len(hull.sample) |
|
1096 |
|
else: |
|
1097 |
|
return hull._npoints |
|
1098 |
|
|
|
1099 |
|
|
|
1100 |
|
|
|
1101 |
|
@property #č stejně nikdo to nebude spouštět |
|
1102 |
|
def nsimplex(hull): |
|
1103 |
|
return hull.sample.nvar * 2 |
|
1104 |
|
|
|
1105 |
|
|
|
1106 |
|
@property |
|
1107 |
|
def A(hull): |
|
1108 |
|
hull._update() |
|
1109 |
|
#č žádná optimizace, ale stejně nikdo to nebude spouštět |
|
1110 |
|
return np.vstack((-hull.orth_basis, hull.orth_basis)) |
|
1111 |
|
|
|
1112 |
|
@property |
|
1113 |
|
def b(hull): |
|
1114 |
|
hull._update() |
|
1115 |
|
return np.concatenate((hull.bm, -hull.bp)) |
|
1116 |
|
|
|
1117 |
|
@property |
|
1118 |
|
def equations(hull): |
|
1119 |
|
hull._update() |
|
1120 |
|
#č žádná optimizace, ale stejně nikdo to nebude spouštět |
|
1121 |
|
A = np.vstack((-hull.orth_basis, hull.orth_basis)) |
|
1122 |
|
b = np.concatenate((hull.bm, -hull.bp)) |
|
1123 |
|
return np.hstack((A,b[:,None])) |
|
1124 |
|
|
|
1125 |
|
def is_inside(hull, nodes): |
|
1126 |
|
hull._update() |
|
1127 |
|
x = nodes.G |
|
1128 |
|
|
|
1129 |
|
#E [normal, offset] forming the hyperplane equation of the facet (see Qhull documentation for more) |
|
1130 |
|
bp = np.atleast_2d(hull.bp).T |
|
1131 |
|
bm = np.atleast_2d(hull.bm).T |
|
1132 |
|
|
|
1133 |
|
# N=ns, E - number of hyperplane equations |
|
1134 |
|
ExN = hull.orth_basis @ x.T |
|
1135 |
|
higher = np.all(ExN < bp, axis=0) |
|
1136 |
|
lower = np.all(ExN > bm, axis=0) |
|
1137 |
|
return np.all((higher, lower), axis=0) |
|
1138 |
|
|
|
1139 |
|
def is_outside(hull, nodes): |
|
1140 |
|
return ~hull.is_inside(nodes) |
|
1141 |
|
|
|
1142 |
|
def get_design_points(hull): |
|
1143 |
|
hull._update() |
|
1144 |
|
sample_model = -hull.A * hull.b.reshape(-1,1) |
|
1145 |
|
return hull.sample.f_model.new_sample(sample_model, space=hull.space) |
|
1146 |
|
|
|
1147 |
|
def get_r(hull): |
|
1148 |
|
hull._update() |
|
1149 |
|
"calculates minimal signed distance to planes. Can therefore be negative" |
|
1150 |
|
return min(-hull.bm[0], hull.bp[0]) |
|
1151 |
|
|
|
1152 |
|
|
|
1153 |
|
|
|
1154 |
|
|
921 |
1155 |
|
|