File convex_hull.py added (mode: 100644) (index 0000000..f1fbb1d) |
|
1 |
|
#!/usr/bin/env python |
|
2 |
|
# coding: utf-8 |
|
3 |
|
|
|
4 |
|
import numpy as np |
|
5 |
|
#from . import IS_stat |
|
6 |
|
from . import sball |
|
7 |
|
#from . import f_models |
|
8 |
|
from scipy import stats |
|
9 |
|
from .candybox import CandyBox |
|
10 |
|
|
|
11 |
|
|
|
12 |
|
#č jako sůl potřebujem statefull třidu, |
|
13 |
|
#č která by schovávala vnitřní implementaciju, |
|
14 |
|
#č která by nabízela jednotné rozhraní pro ostatní kód WellMet |
|
15 |
|
#č a byla by přiměřeně kompatibilní s ConvexHull ze scipy.spatial |
|
16 |
|
#č Ze scipy bych viděl podporu atributů .points, .npoints, .equations |
|
17 |
|
#č Určitě bych chtěl funkce .update() a .is_outside(sample) |
|
18 |
|
#č Atribute .space by ukazoval, v jakém prostoru konvexní obál je vytvořen |
|
19 |
|
|
|
20 |
|
|
|
21 |
|
#ё рельса |
|
22 |
|
#č nepodařílo se mi nějak rozumně zobecnit pro libovolný prostor |
|
23 |
|
#č takže Gauss |
|
24 |
|
#č (každá třida implementuje zvlášť) |
|
25 |
|
#def fire(design_point_G, ns): |
|
26 |
|
# to_fire_G = design_point.G[0] |
|
27 |
|
# r = np.sqrt(np.sum(np.square(to_fire_G))) |
|
28 |
|
# a = to_fire_G / r |
|
29 |
|
# fire_from = stats.norm.cdf(r) |
|
30 |
|
# t = np.linspace(fire_from, 1, ns) |
|
31 |
|
# t = stats.norm.ppf(t) |
|
32 |
|
# fire_G = t.reshape(-1,1) @ a.reshape(1,-1) |
|
33 |
|
# return design_point. |
|
34 |
|
|
|
35 |
|
#č jistě musíme mít nějaký zbytečný kus kódu |
|
36 |
|
#č třida jen pro formu, jen tak na hračku |
|
37 |
|
#č když tečíčky jsou dále jak nejvzdálenější vzorek (bod), |
|
38 |
|
#č tak řekneme, že jsou totálně mimo! |
|
39 |
|
class GBall: |
|
40 |
|
def __init__(hull, sample): |
|
41 |
|
hull.sample = sample |
|
42 |
|
hull.space = 'G' #оӵ малы? Кинлы со кулэ? |
|
43 |
|
|
|
44 |
|
@property |
|
45 |
|
def points(hull): |
|
46 |
|
return hull.sample.G |
|
47 |
|
|
|
48 |
|
@property |
|
49 |
|
def npoints(hull): |
|
50 |
|
return len(hull.sample) |
|
51 |
|
|
|
52 |
|
#def update(hull): pass |
|
53 |
|
|
|
54 |
|
def is_inside(hull, nodes): |
|
55 |
|
R2_hull = np.nanmax(np.sum(np.square(hull.sample.G), axis=1)) |
|
56 |
|
|
|
57 |
|
R2_nodes = np.sum(np.square(nodes.G), axis=1) |
|
58 |
|
return R2_nodes < R2_hull |
|
59 |
|
|
|
60 |
|
def is_outside(hull, nodes): |
|
61 |
|
return ~hull.is_inside(nodes) |
|
62 |
|
|
|
63 |
|
def get_R(hull): |
|
64 |
|
return np.sqrt(np.nanmax(np.sum(np.square(hull.sample.G), axis=1))) |
|
65 |
|
|
|
66 |
|
def get_r(hull): |
|
67 |
|
# zero at best, |
|
68 |
|
# on the safe side, between -R and 0 |
|
69 |
|
#r2_hull = np.nanmin(np.sum(np.square(hull.sample.G), axis=1)) |
|
70 |
|
#return -np.sqrt(r2_hull) |
|
71 |
|
|
|
72 |
|
#č když na get_r nahlížím z hlediska toho, |
|
73 |
|
#č že r vymezuje mezikruží, kde nejsme jistí |
|
74 |
|
#č inside-outside, |
|
75 |
|
#č tak zde get_r vždy musí vracet prostě R |
|
76 |
|
return hull.get_R() |
|
77 |
|
|
|
78 |
|
|
|
79 |
|
def fire(hull, ns): # rail |
|
80 |
|
R2_hull = np.sum(np.square(hull.sample.G), axis=1) |
|
81 |
|
design_index = np.nanargmax(R2_hull) |
|
82 |
|
r = np.sqrt(np.nanmax(R2_hull)) |
|
83 |
|
a = hull.sample.G[design_index] / r |
|
84 |
|
fire_from = stats.norm.cdf(r) |
|
85 |
|
t = np.linspace(fire_from, 1, ns) |
|
86 |
|
t = stats.norm.ppf(t) |
|
87 |
|
fire_G = t.reshape(-1,1) @ a.reshape(1,-1) |
|
88 |
|
|
|
89 |
|
return hull.sample.f_model.new_sample(fire_G, space='G') |
|
90 |
|
|
|
91 |
|
|
|
92 |
|
|
|
93 |
|
#č vyhejbám slovu Box, ať nevnáším ještě víc zmatku v označení |
|
94 |
|
class BrickHull: #č nebo BoundingBrick |
|
95 |
|
def __init__(hull, sample, space='G'): |
|
96 |
|
hull.sample = sample |
|
97 |
|
hull.space = space |
|
98 |
|
|
|
99 |
|
|
|
100 |
|
hull._npoints = 0 |
|
101 |
|
hull.mins = np.full(hull.sample.nvar, np.inf) |
|
102 |
|
hull.maxs = np.full(hull.sample.nvar, -np.inf) |
|
103 |
|
|
|
104 |
|
def _update(hull): |
|
105 |
|
if hull._npoints < hull.npoints: |
|
106 |
|
hull.mins = np.nanmin(hull.points, axis=0) |
|
107 |
|
hull.maxs = np.nanmax(hull.points, axis=0) |
|
108 |
|
hull._npoints = hull.npoints |
|
109 |
|
|
|
110 |
|
|
|
111 |
|
def is_inside(hull, nodes): |
|
112 |
|
hull._update() |
|
113 |
|
x = getattr(nodes, hull.space) |
|
114 |
|
more = np.all(x < hull.maxs, axis=1) |
|
115 |
|
less = np.all(x > hull.mins, axis=1) |
|
116 |
|
|
|
117 |
|
return np.all((more, less), axis=0) #np.all(np.array((more, less)), axis=0) |
|
118 |
|
|
|
119 |
|
def is_outside(hull, nodes): |
|
120 |
|
return ~hull.is_inside(nodes) |
|
121 |
|
|
|
122 |
|
def get_design_points(hull): |
|
123 |
|
hull._update() |
|
124 |
|
#sample_model = -hull.A * hull.b.reshape(-1,1) |
|
125 |
|
sample_model = np.vstack((np.diag(hull.maxs), np.diag(hull.mins))) |
|
126 |
|
return hull.sample.f_model.new_sample(sample_model, space=hull.space) |
|
127 |
|
|
|
128 |
|
# shortcut for Ghull |
|
129 |
|
# valid only if space==G |
|
130 |
|
def get_r(hull): |
|
131 |
|
if hull.space=='G': |
|
132 |
|
hull._update() |
|
133 |
|
return np.min((-hull.mins, hull.maxs)) |
|
134 |
|
else: |
|
135 |
|
return 0 |
|
136 |
|
|
|
137 |
|
@property |
|
138 |
|
def points(hull): |
|
139 |
|
return getattr(hull.sample, hull.space) |
|
140 |
|
|
|
141 |
|
@property |
|
142 |
|
def npoints(hull): |
|
143 |
|
return len(hull.sample) |
|
144 |
|
|
|
145 |
|
@property |
|
146 |
|
def A(hull): |
|
147 |
|
hull._update() |
|
148 |
|
#č žádná optimizace, ale kdo tu funkci bude spouštět? |
|
149 |
|
diag = np.diag(np.ones(hull.sample.nvar)) |
|
150 |
|
return np.vstack((diag, -diag)) |
|
151 |
|
|
|
152 |
|
@property |
|
153 |
|
def b(hull): |
|
154 |
|
hull._update() |
|
155 |
|
return np.concatenate((hull.maxs, -hull.mins)) |
|
156 |
|
|
|
157 |
|
@property |
|
158 |
|
def equations(hull): |
|
159 |
|
hull._update() |
|
160 |
|
#č žádná optimizace, ale kdo tu funkci bude spouštět? |
|
161 |
|
diag = np.diag(np.ones(hull.sample.nvar)) |
|
162 |
|
|
|
163 |
|
A = np.vstack((diag, -diag)) |
|
164 |
|
b = np.concatenate((-hull.maxs, hull.mins)) |
|
165 |
|
return np.hstack((A,b[:,None])) |
|
166 |
|
|
|
167 |
|
def fire(hull, ns): # boom |
|
168 |
|
sample_U = hull.sample.U |
|
169 |
|
|
|
170 |
|
U_mins = np.nanmin(sample_U, axis=0) |
|
171 |
|
U_maxs = np.nanmax(sample_U, axis=0) |
|
172 |
|
|
|
173 |
|
#č min nebo max? |
|
174 |
|
if np.nanmax(U_mins) > (1-np.nanmin(U_maxs)): |
|
175 |
|
#č expandujeme se dolu |
|
176 |
|
var = np.nanargmax(U_mins) |
|
177 |
|
value = U_mins[var] |
|
178 |
|
t = np.linspace(value, 0, ns, endpoint=False) |
|
179 |
|
else: |
|
180 |
|
#č expandujeme se nahoru |
|
181 |
|
var = np.nanargmin(U_maxs) |
|
182 |
|
value = U_maxs[var] |
|
183 |
|
t = np.linspace(value, 1, ns, endpoint=False) |
|
184 |
|
|
|
185 |
|
boom = np.random.random((ns, hull.sample.nvar)) |
|
186 |
|
boom[:, var] = t |
|
187 |
|
|
|
188 |
|
return hull.sample.f_model.new_sample(boom, space='U') |
|
189 |
|
|
|
190 |
|
|
|
191 |
|
class DirectHull: |
|
192 |
|
def __init__(hull, sample, direct_plan, space='G'): |
|
193 |
|
hull.sample = sample |
|
194 |
|
hull.direct_plan = direct_plan |
|
195 |
|
hull.space = space |
|
196 |
|
|
|
197 |
|
hull.regen() |
|
198 |
|
|
|
199 |
|
def regen(hull): |
|
200 |
|
hull._npoints = 0 |
|
201 |
|
hull.mins = np.full(hull.sample.nvar, np.inf) |
|
202 |
|
hull.maxs = np.full(hull.sample.nvar, -np.inf) |
|
203 |
|
|
|
204 |
|
hull.bp = np.full(len(hull.direct_plan), -np.inf) |
|
205 |
|
hull.bm = np.full(len(hull.direct_plan), np.inf) |
|
206 |
|
|
|
207 |
|
#hull.update() |
|
208 |
|
|
|
209 |
|
#č nejsem jist, jestli ten update vůbec dělat. |
|
210 |
|
#č lze navrhnout třidu aj tak, že sama bude hlídat změny. |
|
211 |
|
#č Jenomže, co kdybychom ten automatický update nechtěli? |
|
212 |
|
def _update(hull): |
|
213 |
|
if hull._npoints < hull.npoints: |
|
214 |
|
#hull.points = getattr(hull.sample, hull.space) |
|
215 |
|
new_points = hull.points[hull._npoints:] |
|
216 |
|
|
|
217 |
|
new_mins = np.nanmin(new_points, axis=0) |
|
218 |
|
new_maxs = np.nanmax(new_points, axis=0) |
|
219 |
|
|
|
220 |
|
hull.mins = np.nanmin((new_mins, hull.mins), axis=0) |
|
221 |
|
hull.maxs = np.nanmax((new_maxs, hull.maxs), axis=0) |
|
222 |
|
|
|
223 |
|
A = hull.direct_plan @ new_points.T |
|
224 |
|
new_bp = np.nanmax(A, axis=1) |
|
225 |
|
new_bm = np.nanmin(A, axis=1) |
|
226 |
|
|
|
227 |
|
hull.bp = np.nanmax((new_bp, hull.bp), axis=0) |
|
228 |
|
hull.bm = np.nanmin((new_bm, hull.bm), axis=0) |
|
229 |
|
|
|
230 |
|
hull._npoints = hull.npoints |
|
231 |
|
|
|
232 |
|
@property |
|
233 |
|
def points(hull): |
|
234 |
|
return getattr(hull.sample, hull.space) |
|
235 |
|
|
|
236 |
|
@property |
|
237 |
|
def npoints(hull): |
|
238 |
|
return len(hull.sample) |
|
239 |
|
|
|
240 |
|
@property |
|
241 |
|
def A(hull): |
|
242 |
|
hull._update() |
|
243 |
|
#č žádná optimizace, ale kdo tu funkci bude spouštět? |
|
244 |
|
diag = np.diag(np.ones(hull.sample.nvar)) |
|
245 |
|
return np.vstack((diag, -diag, -hull.direct_plan, hull.direct_plan)) |
|
246 |
|
|
|
247 |
|
@property |
|
248 |
|
def b(hull): |
|
249 |
|
hull._update() |
|
250 |
|
return np.concatenate((-hull.maxs, hull.mins, hull.bm, -hull.bp)) |
|
251 |
|
|
|
252 |
|
@property |
|
253 |
|
def equations(hull): |
|
254 |
|
hull._update() |
|
255 |
|
#č žádná optimizace, ale kdo tu funkci bude spouštět? |
|
256 |
|
diag = np.diag(np.ones(hull.sample.nvar)) |
|
257 |
|
|
|
258 |
|
A = np.vstack((diag, -diag, -hull.direct_plan, hull.direct_plan)) |
|
259 |
|
b = np.concatenate((-hull.maxs, hull.mins, hull.bm, -hull.bp)) |
|
260 |
|
return np.hstack((A,b[:,None])) |
|
261 |
|
|
|
262 |
|
def is_inside(hull, nodes): |
|
263 |
|
hull._update() |
|
264 |
|
x = getattr(nodes, hull.space) |
|
265 |
|
|
|
266 |
|
more = np.all(x < hull.maxs, axis=1) |
|
267 |
|
less = np.all(x > hull.mins, axis=1) |
|
268 |
|
|
|
269 |
|
|
|
270 |
|
#E [normal, offset] forming the hyperplane equation of the facet (see Qhull documentation for more) |
|
271 |
|
A = hull.direct_plan |
|
272 |
|
bp = np.atleast_2d(hull.bp).T |
|
273 |
|
bm = np.atleast_2d(hull.bm).T |
|
274 |
|
|
|
275 |
|
# N=ns, E - number of hyperplane equations |
|
276 |
|
ExN = A @ x.T |
|
277 |
|
higher = np.all(ExN < bp, axis=0) |
|
278 |
|
lower = np.all(ExN > bm, axis=0) |
|
279 |
|
return np.all((more, less, higher, lower), axis=0) |
|
280 |
|
|
|
281 |
|
def is_outside(hull, nodes): |
|
282 |
|
return ~hull.is_inside(nodes) |
|
283 |
|
|
|
284 |
|
def get_design_points(hull): |
|
285 |
|
hull._update() |
|
286 |
|
sample_model = -hull.A * hull.b.reshape(-1,1) |
|
287 |
|
return hull.sample.f_model.new_sample(sample_model, space=hull.space) |
|
288 |
|
|
|
289 |
|
# shortcut for Ghull |
|
290 |
|
# valid only if space==G |
|
291 |
|
def get_r(hull): |
|
292 |
|
if hull.space=='G': |
|
293 |
|
hull._update() |
|
294 |
|
return min(-np.max(hull.mins), np.min(hull.maxs), -np.max(hull.bm), np.min(hull.bp)) |
|
295 |
|
else: |
|
296 |
|
return 0 |
|
297 |
|
|
|
298 |
|
def fire(hull, ns): |
|
299 |
|
if hull.space == 'G': |
|
300 |
|
A = hull.equations[:,:-1] |
|
301 |
|
b = hull.equations[:,-1] |
|
302 |
|
|
|
303 |
|
to_fire = np.argmax(b) |
|
304 |
|
a = A[to_fire] |
|
305 |
|
fire_from = stats.norm.cdf(hull.get_r()) |
|
306 |
|
t = np.linspace(fire_from, 1, ns) |
|
307 |
|
t = stats.norm.ppf(t) |
|
308 |
|
fire_G = t.reshape(-1,1) @ a.reshape(1,-1) |
|
309 |
|
|
|
310 |
|
return hull.sample.f_model.new_sample(fire_G, space='G') |
|
311 |
|
|
|
312 |
|
#č mým úkolem při návrhu této třidy je pořádně všecko zkomplikovat. |
|
313 |
|
#č Dostávám za to peníze. |
|
314 |
|
#č Takže. Udělám 3 druhů estimátorů |
|
315 |
|
# convex_hull_estimation -2: inside, -1: outside |
|
316 |
|
# shell_estimation -22: inner, -3: shell, -11: outer |
|
317 |
|
# ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer |
|
318 |
|
class Ghull: |
|
319 |
|
def __init__(self, hull, design=None): |
|
320 |
|
self.hull = hull |
|
321 |
|
self.design = design |
|
322 |
|
self.shell = sball.Shell(hull.sample.nvar) |
|
323 |
|
self.sample = hull.sample |
|
324 |
|
|
|
325 |
|
# explicit wrapping |
|
326 |
|
self.fire = hull.fire |
|
327 |
|
|
|
328 |
|
def get_R(self): |
|
329 |
|
sum_squared = np.sum(np.square(self.sample.G), axis=1) |
|
330 |
|
#index = np.argmax(sum_squared) |
|
331 |
|
return np.sqrt(np.nanmax(sum_squared)) |
|
332 |
|
|
|
333 |
|
def get_r(self): |
|
334 |
|
"calculates minimal signed distance to planes. Can therefore be negative" |
|
335 |
|
#return -np.nanmax(self.hull.b) |
|
336 |
|
return self.hull.get_r() |
|
337 |
|
|
|
338 |
|
def get_shell_estimation(self): |
|
339 |
|
shell = self.shell |
|
340 |
|
r = self.get_r() |
|
341 |
|
R = self.get_R() |
|
342 |
|
|
|
343 |
|
if r<0: |
|
344 |
|
shell.set_bounds(0, R) |
|
345 |
|
else: |
|
346 |
|
shell.set_bounds(r, R) |
|
347 |
|
|
|
348 |
|
# shell_estimation -22: inner, -3: shell, -11: outer |
|
349 |
|
shell_estimation = {-22:shell.ps, -3: shell.p_shell, -11: shell.pf} |
|
350 |
|
global_stats = {"r":r, "R":R, "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf} |
|
351 |
|
return shell_estimation, global_stats |
|
352 |
|
|
|
353 |
|
def integrate(self, nis): |
|
354 |
|
#č no to teda disajn třidy je doopravdy hroznej |
|
355 |
|
# .get_shell_estimation() will calculate radiuses and will update shell |
|
356 |
|
shell_estimation, global_stats = self.get_shell_estimation() |
|
357 |
|
nodes_G = self.shell.rvs(nis) |
|
358 |
|
nodes = self.sample.f_model.new_sample(nodes_G, space='G') |
|
359 |
|
mask = self.hull.is_outside(nodes) |
|
360 |
|
#č nevím proč, ale kdysi mě to vyšlo rychlejc jak obyčejný součet |
|
361 |
|
nf = len(mask[mask]) |
|
362 |
|
ns = nis - nf |
|
363 |
|
# shell_estimation -22: inner, -3: shell, -11: outer |
|
364 |
|
p_shell = shell_estimation.pop(-3) |
|
365 |
|
shell_pf = nf/nis * p_shell |
|
366 |
|
shell_ps = ns/nis * p_shell |
|
367 |
|
|
|
368 |
|
# ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer |
|
369 |
|
ghull_estimation = shell_estimation; del(shell_estimation) |
|
370 |
|
ghull_estimation[-21] = shell_ps |
|
371 |
|
ghull_estimation[-12] = shell_pf |
|
372 |
|
global_stats["shell inside"] = shell_ps |
|
373 |
|
global_stats["shell outside"] = shell_pf |
|
374 |
|
|
|
375 |
|
# convex_hull_estimation -2: inside, -1: outside |
|
376 |
|
inside = shell_ps + self.shell.ps |
|
377 |
|
outside = shell_pf + self.shell.pf |
|
378 |
|
convex_hull_estimation = {-2: inside, -1: outside} |
|
379 |
|
global_stats["inside"] = inside |
|
380 |
|
global_stats["outside"] = outside |
|
381 |
|
|
|
382 |
|
#nodes = self.sample.f_model.new_sample(nodes_G, space='G') |
|
383 |
|
# -2 = 'inside' -1 = 'outside' |
|
384 |
|
nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask) |
|
385 |
|
|
|
386 |
|
return nodes, ghull_estimation, convex_hull_estimation, global_stats |
|
387 |
|
|
|
388 |
|
|
|
389 |
|
|
|
390 |
|
|
|
391 |
|
##č napadlo mě zababáchnuť třidu, která by se sama starala o všem co se tyče |
|
392 |
|
##č vnější domény. Nešlo mě totíž to udělat jednou funkcí, bylo by velmi |
|
393 |
|
##č špatné z hlediska zodpovednosti kódu. Tak to všecko zabalíme to třidy |
|
394 |
|
##č a odebereme z už beztak přetíženého blackboxu část komplexity |
|
395 |
|
## keywords: ISSI, estimation, outside, ConvexHull, Sball, IS kolem středních hodnot |
|
396 |
|
#class Shull: # issi_estimate_outside |
|
397 |
|
# def __init__(sx, f, model_space, sampling_space='G', \ |
|
398 |
|
# powerset_correction=True, incremental=True, design=None): |
|
399 |
|
# #č tím powerset_corretion je myšlena úplná soustava jevů, |
|
400 |
|
# #č tj. vyrovnaní s použitím míry vnější i vnitřní |
|
401 |
|
# #č powerset_correction=True přídá -2 (inside) jev do ISSI |
|
402 |
|
# |
|
403 |
|
# #č zde f-ko musí taky obsahovat vzorky! |
|
404 |
|
# sx.f = f |
|
405 |
|
# sx.model_space = model_space |
|
406 |
|
# |
|
407 |
|
# #č kašlu na to, pokud uživatel zadal nesmysl, tak |
|
408 |
|
# #č sam převedu do nečěho smyslúplnějšího |
|
409 |
|
# _dict = {'R':'Rn', 'aR':'aRn', 'P':'GK', 'aP':'aGK', 'U':'G', 'aU':'aG'} |
|
410 |
|
# if sampling_space in _dict: |
|
411 |
|
# sx.sampling_space = _dict[sampling_space] |
|
412 |
|
# else: |
|
413 |
|
# sx.sampling_space = sampling_space |
|
414 |
|
# |
|
415 |
|
# sx.design = design |
|
416 |
|
# |
|
417 |
|
# sampled_plan_model = getattr(f, model_space) |
|
418 |
|
# #č žádná kontrola chyb - nechť to spadné, když má spadnout! |
|
419 |
|
# sx.convex_hull = spatial.ConvexHull(sampled_plan_model, incremental=incremental) |
|
420 |
|
# |
|
421 |
|
# # current outside probability estimation |
|
422 |
|
# sx.p_out = 0.5 # implicit value |
|
423 |
|
# sx.sball = sball.Sball(f.nvar) |
|
424 |
|
# sx.base_r = sx.sball.get_r(0.5) # we want in average 50/50 ratio |
|
425 |
|
# |
|
426 |
|
# |
|
427 |
|
# # -1 = 'outside' |
|
428 |
|
# # -2 = 'inside' |
|
429 |
|
# #č založme ISSI |
|
430 |
|
# sx.powerset_correction = powerset_correction |
|
431 |
|
# #č potřebuji pro korektnost mít před integrací zadané jevy |
|
432 |
|
# if powerset_correction: |
|
433 |
|
# sx.oiss = IS_stat.ISSI([-1, -2]) |
|
434 |
|
# else: |
|
435 |
|
# sx.oiss = IS_stat.ISSI([-1]) |
|
436 |
|
# |
|
437 |
|
# |
|
438 |
|
# |
|
439 |
|
# def increment(sx, input_sample): |
|
440 |
|
# #č sample by měl byt jíž převeden na f (v .add_sample()), |
|
441 |
|
# #č zodpovidá za to volajicí kód! |
|
442 |
|
# sx.convex_hull.add_points(getattr(input_sample, sx.model_space)) |
|
443 |
|
# |
|
444 |
|
# |
|
445 |
|
# |
|
446 |
|
# # require |
|
447 |
|
# def integrate(sx, nis): |
|
448 |
|
# # getting rid of the old estimations |
|
449 |
|
# sx.oiss.delete_event_data(-1) |
|
450 |
|
# sx.oiss.events.append(-1) #č už tak trošku sahám do vnitřku cizí třidy |
|
451 |
|
# if sx.powerset_correction: |
|
452 |
|
# sx.oiss.delete_event_data(-2) |
|
453 |
|
# sx.oiss.events.append(-2) #č a záse |
|
454 |
|
# |
|
455 |
|
# #č posunutí středu trošiňku porušuje předpoklady, za kterých |
|
456 |
|
# #č Sball volí rozdělení, ale přečo se mi stavá, |
|
457 |
|
# #č že ve více dimenzích Shull na začatku prostě nemůže |
|
458 |
|
# #č trefit ten blbej... v podstatě simplex |
|
459 |
|
# #č Těžiště, přesnějí, prostě nějaký střed můžeme najít dvěma způsoby: |
|
460 |
|
# #č 1. mean(vertices.model_space).sampling_space |
|
461 |
|
# #č 2. mean(vertices.sampling_space) |
|
462 |
|
# #č myslím si, že ten první (a bez váh) je stabilnější |
|
463 |
|
# #č (víme, jak střed vrcholů v nějakém prostoru může ani netrefit do |
|
464 |
|
# #č simplexu v původním prostoru) |
|
465 |
|
## vertices_model = sx.convex_hull.points[sx.convex_hull.vertices] |
|
466 |
|
## barycenter_model = np.mean(vertices_model, axis=0) |
|
467 |
|
## if sx.model_space == sx.sampling_space: |
|
468 |
|
## sx.barycenter_sampling = barycenter_model |
|
469 |
|
## else: |
|
470 |
|
## barycenter = sx.f.new_sample(barycenter_model, sx.model_space) |
|
471 |
|
## sx.barycenter_sampling = np.squeeze(getattr(barycenter, sx.sampling_space)) |
|
472 |
|
# |
|
473 |
|
# #č rozhodl jsem, že shull musí odhadovat outside, ne motat se centrem Brna. |
|
474 |
|
# #č předpokladám, že uživatel zadal buď G, nebo aspoň Rn prostor |
|
475 |
|
# #sx.barycenter_sampling = np.full(sx.f.nvar, 0, dtype=np.int8) |
|
476 |
|
# |
|
477 |
|
# # first step |
|
478 |
|
# nodes = sx._sample_sball(nis) |
|
479 |
|
# mask = nodes.is_outside |
|
480 |
|
# |
|
481 |
|
# |
|
482 |
|
# cut_off_out = int(nis/3) |
|
483 |
|
# cut_off_in = int(nis/3) |
|
484 |
|
# #č robím cyklus dokud nesberu dostatečně teček. |
|
485 |
|
# #č To je fakt nejrobustnější řešení, co mě napadá |
|
486 |
|
# # while (number_of_out_nodes or number_of_nodes_inside is too_small) |
|
487 |
|
# while (len(mask[mask]) < cut_off_out): # or (len(mask[~mask]) < cut_off_in): |
|
488 |
|
# print(sx.__class__.__name__ + ":", "cut off! Outsides: %s, insides: %s, p_out=%s"\ |
|
489 |
|
# % (len(mask[mask]), len(mask[~mask]), sx.p_out)) |
|
490 |
|
# #č je třeba jenom sehnat dostatečně bodíků a utikat |
|
491 |
|
# nodes.add_sample(sx._sample_sball(nis)) |
|
492 |
|
# mask = nodes.is_outside |
|
493 |
|
# |
|
494 |
|
# #č když neprovadíme výrovnání, tak ten vnitršek nachren nepotřebujem |
|
495 |
|
# if (len(mask[~mask]) < cut_off_in) and sx.powerset_correction: |
|
496 |
|
# print(sx.__class__.__name__ + ":", \ |
|
497 |
|
# "cut off inside (%s of %s needed)! Do simplex-based integration of convex hull"\ |
|
498 |
|
# % (len(mask[~mask]), cut_off_in)) |
|
499 |
|
# |
|
500 |
|
# nodes.add_sample(sx._sample_simplex(nis)) |
|
501 |
|
# |
|
502 |
|
# |
|
503 |
|
# return nodes |
|
504 |
|
# |
|
505 |
|
# |
|
506 |
|
# |
|
507 |
|
# def _sample_simplex(sx, nis): |
|
508 |
|
# |
|
509 |
|
# #č f-ko sice musí odkazovat na aktuální f_model |
|
510 |
|
# #č ale na druhou stranu normálně ._integrate_simplex() |
|
511 |
|
# #č potřebujeme pouze jednou - hned při vytvaření |
|
512 |
|
# vertices = sx.f[sx.convex_hull.vertices] |
|
513 |
|
# nvar = vertices.nvar |
|
514 |
|
# |
|
515 |
|
# # IS_like uses .new_sample method, so vertices can not be a SampleBox object |
|
516 |
|
# # |
|
517 |
|
# #č IS_like těžiště počítá v sampling_space, ale asi mu to až tak nevadí |
|
518 |
|
# # |
|
519 |
|
# # already divided by nsim in variance formule |
|
520 |
|
# # divide by /(nvar+1)/(nvar+2) from simplex inertia tensor solution |
|
521 |
|
# # multiply by simplex_volume, but it looks like it shouldn't be here |
|
522 |
|
# # for simplex: d = nvar+2 |
|
523 |
|
# #č sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku |
|
524 |
|
# nodes = IS_stat.IS_like(vertices, sampling_space=sx.sampling_space, \ |
|
525 |
|
# nis=nis, d=nvar+2, design=sx.design) |
|
526 |
|
# |
|
527 |
|
# #č indikatorová funkce |
|
528 |
|
# sx.is_outside(nodes) |
|
529 |
|
# |
|
530 |
|
# # for IS_stats |
|
531 |
|
# #č zkoušel jsem zadavat celou serii - zhoršuje to odhady |
|
532 |
|
# #č nemůžeme důvěrovat tomu, jak ten slepej simplex vidí svět |
|
533 |
|
# weights = nodes.w[~nodes.is_outside] |
|
534 |
|
# sx.oiss.add_single_event_data(weights, event=-2, nis=nis) |
|
535 |
|
# # add_single_event_data() do not calculate estimations itself |
|
536 |
|
# sx.oiss.get_estimations() |
|
537 |
|
# return nodes |
|
538 |
|
# |
|
539 |
|
# |
|
540 |
|
# |
|
541 |
|
# def _sample_sball(sx, nis): |
|
542 |
|
# nvar = sx.f.nvar |
|
543 |
|
# sampling_r, sx.p_out = sx.sball.get_r_iteration(sx.p_out) |
|
544 |
|
# #sampling_r = sx.sball.get_r(sx.p_out) |
|
545 |
|
# |
|
546 |
|
# #č asi tam bylo sampling_r/base_r, že? |
|
547 |
|
# std_ball = sampling_r/sx.base_r |
|
548 |
|
# #č chcu std=1, když p_out -> 1 |
|
549 |
|
# #č a std=sball_solušn, když p_out -> 0 |
|
550 |
|
# #č surovější nevymyslíš! |
|
551 |
|
# std = std_ball + sx.p_out |
|
552 |
|
# #č u stats.norm zadáváme směrodatnou odchylku, je to asi správné |
|
553 |
|
# #h = f_models.UnCorD([stats.norm(sx.barycenter_sampling[i], std) for i in range(nvar)]) |
|
554 |
|
# #nodes = IS_stat.IS(sx.f, h, space_from_h='R', space_to_f=sx.sampling_space, Nsim=nis) |
|
555 |
|
# |
|
556 |
|
# #norm_params = [(sx.barycenter_sampling[i], std) for i in range(nvar)] |
|
557 |
|
# nodes = IS_stat.IS_norm(sx.f, std=std, \ |
|
558 |
|
# sampling_space=sx.sampling_space, nis=nis, design=sx.design) |
|
559 |
|
# |
|
560 |
|
# outside_measure = sx._apply_nodes(nodes, nis) |
|
561 |
|
# |
|
562 |
|
# #č pro přiště |
|
563 |
|
# sx.p_out = (sx.p_out + outside_measure) / 2 |
|
564 |
|
# |
|
565 |
|
# return nodes |
|
566 |
|
# |
|
567 |
|
# |
|
568 |
|
# def _apply_nodes(sx, nodes, nis): |
|
569 |
|
# #č indikatorová funkce |
|
570 |
|
# sx.is_outside(nodes) |
|
571 |
|
# |
|
572 |
|
# # for IS_stats |
|
573 |
|
# if sx.powerset_correction: |
|
574 |
|
# #č získáme výrovnaný odhad - je to skoro zdarma |
|
575 |
|
# #svar = (sampling_r/sx.base_r)**2 # svar like sampling_variance |
|
576 |
|
# #č kdysi snažil jsem něco odvést, moc se mi to nepovedlo |
|
577 |
|
# #č je to jen tak, jeden z pokusu, hrubej nastřel |
|
578 |
|
# #im = svar**nvar * np.exp(nvar/svar - nvar) |
|
579 |
|
# |
|
580 |
|
# #č radší ne. IM špatně zachycuje nizkou důvěru k tomu, co nemá vlastní tečky |
|
581 |
|
# sx.oiss.add_IS_serie(nodes.w, nodes.event_id, implicit_multiplicator=np.inf) |
|
582 |
|
# outside_measure = sx.oiss.estimations[-1] |
|
583 |
|
# else: |
|
584 |
|
# weights = nodes.w[nodes.is_outside] |
|
585 |
|
# #č IM všecko pokazí, jakmile začnu přídávát další jevy |
|
586 |
|
# sx.oiss.add_single_event_data(weights, event=-1, nis=nis) |
|
587 |
|
# # add_single_event_data() do not calculate estimations itself |
|
588 |
|
# weighted_mean, __, events = sx.oiss.get_means() |
|
589 |
|
# # outside probability |
|
590 |
|
# #č nevyrovnané! |
|
591 |
|
# outside_measure = weighted_mean[events==-1][0] |
|
592 |
|
# |
|
593 |
|
# return outside_measure |
|
594 |
|
# |
|
595 |
|
# |
|
596 |
|
# |
|
597 |
|
# def is_outside(sx, nodes): |
|
598 |
|
# node_coordinates = getattr(nodes, sx.model_space) |
|
599 |
|
# mask = is_outside(sx.convex_hull, node_coordinates) |
|
600 |
|
# |
|
601 |
|
# # -1 = 'outside' |
|
602 |
|
# # -2 = 'inside' |
|
603 |
|
# nodes.event_id = mask - 2 |
|
604 |
|
# nodes.is_outside = mask |
|
605 |
|
# |
|
606 |
|
# |
|
607 |
|
|
|
608 |
|
|
File qt_plot.py changed (mode: 100644) (index d3bfdde..027fd3f) |
... |
... |
from pyqtgraph import console |
10 |
10 |
import numpy as np |
import numpy as np |
11 |
11 |
import pandas as pd # required for estimation graph |
import pandas as pd # required for estimation graph |
12 |
12 |
|
|
|
13 |
|
#ё я сдаюсь. |
|
14 |
|
#č quadpy tak se stavá povinnou závislostí |
|
15 |
|
#č potrebuju pro HullEstimation widget |
|
16 |
|
import quadpy |
|
17 |
|
|
|
18 |
|
|
13 |
19 |
from . import estimation as stm |
from . import estimation as stm |
14 |
20 |
from . import misc |
from . import misc |
15 |
21 |
from . import stm_df |
from . import stm_df |
|
22 |
|
from . import sball |
|
23 |
|
from . import schemes |
|
24 |
|
from . import convex_hull as khull |
|
25 |
|
|
16 |
26 |
|
|
17 |
27 |
def qt_gui_plot_2d(sample_box, space='R', *args, **kwargs): |
def qt_gui_plot_2d(sample_box, space='R', *args, **kwargs): |
18 |
28 |
""" |
""" |
|
... |
... |
class QtGuiPlot2D(QtGui.QMainWindow): |
41 |
51 |
super().__init__() |
super().__init__() |
42 |
52 |
self.setWindowTitle("%sD: %s" %(sample_box.nvar, sample_box.gm_signature)) |
self.setWindowTitle("%sD: %s" %(sample_box.nvar, sample_box.gm_signature)) |
43 |
53 |
|
|
|
54 |
|
# for debug |
|
55 |
|
# container for errors |
|
56 |
|
# to trace errors |
|
57 |
|
self.errors = [] |
|
58 |
|
|
44 |
59 |
self.kwargs = kwargs |
self.kwargs = kwargs |
45 |
60 |
self.sample_box = sample_box |
self.sample_box = sample_box |
46 |
61 |
#sample_box.sample_box._log = self.logger |
#sample_box.sample_box._log = self.logger |
|
... |
... |
class Giracles(Series): |
488 |
503 |
f_model = self.w.sample_box.f_model |
f_model = self.w.sample_box.f_model |
489 |
504 |
sample_G = self.prebound * r |
sample_G = self.prebound * r |
490 |
505 |
sample = f_model.new_sample(sample_G, space='G', extend=True) |
sample = f_model.new_sample(sample_G, space='G', extend=True) |
491 |
|
return self.add_serie(sample, index=None, **plot_kwargs) |
|
|
506 |
|
return self.add_serie(sample, index=index, **plot_kwargs) |
492 |
507 |
|
|
493 |
508 |
|
|
|
509 |
|
|
|
510 |
|
|
|
511 |
|
class InfiniteLines(Series): |
|
512 |
|
|
|
513 |
|
def add_line(self, space='G', index=None, **plot_kwargs): |
|
514 |
|
plot_item = self.w.central_widget.addLine(**plot_kwargs) |
|
515 |
|
if space == self.w.space: |
|
516 |
|
plot_item.show() |
|
517 |
|
else: |
|
518 |
|
plot_item.hide() |
|
519 |
|
|
|
520 |
|
if index is None: |
|
521 |
|
index = len(self.items) |
|
522 |
|
elif index in self.items: |
|
523 |
|
# kind of update, then |
|
524 |
|
#č musíme korektně odebrat předchozí kresbu |
|
525 |
|
self.remove_item(index) |
|
526 |
|
|
|
527 |
|
self.items[index] = [space, plot_item, plot_kwargs] |
|
528 |
|
return plot_item |
|
529 |
|
|
|
530 |
|
def _redraw(self): |
|
531 |
|
if self.autoredraw: |
|
532 |
|
for item in self.items.values(): |
|
533 |
|
space, _invalid_plot_item, plot_dict = item |
|
534 |
|
item[1] = self.w.central_widget.addLine(**plot_dict) |
|
535 |
|
else: |
|
536 |
|
self.items.clear() |
|
537 |
|
|
|
538 |
|
|
|
539 |
|
def space_update(self): |
|
540 |
|
for item in self.items.values(): |
|
541 |
|
space, plot_item, __plot_dict = item |
|
542 |
|
if space == self.w.space: |
|
543 |
|
plot_item.show() |
|
544 |
|
else: |
|
545 |
|
plot_item.hide() |
|
546 |
|
|
|
547 |
|
|
494 |
548 |
""" |
""" |
495 |
549 |
============== |
============== |
496 |
550 |
у График люкет |
у График люкет |
|
... |
... |
class HullEstimationWidget(pg.LayoutWidget): |
2804 |
2858 |
self.sb_item = samplebox_item |
self.sb_item = samplebox_item |
2805 |
2859 |
|
|
2806 |
2860 |
self.giracle = Giracles(w=samplebox_item, autoredraw=False, nrod=200) |
self.giracle = Giracles(w=samplebox_item, autoredraw=False, nrod=200) |
|
2861 |
|
#ё hyperplanes? Кого ты обманываешь, Alexi? |
|
2862 |
|
#č tak. už je to equation_planes :) |
|
2863 |
|
self.equation_planes = InfiniteLines(w=samplebox_item, autoredraw=False) |
|
2864 |
|
|
|
2865 |
|
# signals handling: |
|
2866 |
|
# box_runned - handled by .on_box_run() method |
|
2867 |
|
# space_changed - handled automatically by smart items |
|
2868 |
|
# slice_changed - setted up clear() on smart items |
|
2869 |
|
# redraw_called - handled automatically by smart items |
|
2870 |
|
# estimation_added - class does not handle the signal |
|
2871 |
|
# nor emits the signal itself |
|
2872 |
|
# (as does not save estimations anywhere for now) |
|
2873 |
|
|
|
2874 |
|
# todo: design estimation record |
2807 |
2875 |
|
|
2808 |
2876 |
self.sb_item.box_runned.connect(self.on_box_run) |
self.sb_item.box_runned.connect(self.on_box_run) |
2809 |
|
#č Alexi, ten signal může té funkce posilát bůhví co navíc. Uvidíme. |
|
|
2877 |
|
#č Alexi, ten signal může té funkce posilát bůhví co navíc. |
|
2878 |
|
#č (ano. clear() nebere žádné argumenty, takže v cajku) |
2810 |
2879 |
self.sb_item.slice_changed.connect(self.giracle.clear) |
self.sb_item.slice_changed.connect(self.giracle.clear) |
2811 |
|
#self.sb_item.space_changed.connect(self.on_space_changed) |
|
2812 |
|
#self.sb_item.redraw_called.connect(self.redraw) |
|
|
2880 |
|
self.sb_item.slice_changed.connect(self.equation_planes.clear) |
|
2881 |
|
#č ty naše chytré prvky giracle a equation_planes |
|
2882 |
|
#č hlídají space_changed a redraw_called sami. |
2813 |
2883 |
|
|
|
2884 |
|
self.schemes = dict() |
|
2885 |
|
self.ndim = self.sb_item.sample_box.nvar |
2814 |
2886 |
#☺ na internetu všichni tak dělaj |
#☺ na internetu všichni tak dělaj |
2815 |
2887 |
self.setup() |
self.setup() |
2816 |
2888 |
|
|
|
... |
... |
class HullEstimationWidget(pg.LayoutWidget): |
2832 |
2904 |
#size_policy.setHorizontalPolicy(QtGui.QSizePolicy.Expanding) |
#size_policy.setHorizontalPolicy(QtGui.QSizePolicy.Expanding) |
2833 |
2905 |
#self.toolbar.setSizePolicy(size_policy) |
#self.toolbar.setSizePolicy(size_policy) |
2834 |
2906 |
|
|
|
2907 |
|
#č draw_convex_hull ja navržena tak, aby brala jíž hotový hull |
|
2908 |
|
# if self.ndim == 2: |
|
2909 |
|
# action = self.toolbar.addAction("draw convex hull", self.draw_convex_hull) |
|
2910 |
|
# btn = self.toolbar.widgetForAction(action) |
|
2911 |
|
# btn.setAutoRaise(False) |
|
2912 |
|
|
2835 |
2913 |
action = self.toolbar.addAction("shell out!", self.get_shell_estimation) |
action = self.toolbar.addAction("shell out!", self.get_shell_estimation) |
2836 |
2914 |
btn = self.toolbar.widgetForAction(action) |
btn = self.toolbar.widgetForAction(action) |
2837 |
2915 |
btn.setAutoRaise(False) |
btn.setAutoRaise(False) |
|
... |
... |
class HullEstimationWidget(pg.LayoutWidget): |
2847 |
2925 |
btn.setAutoRaise(False) |
btn.setAutoRaise(False) |
2848 |
2926 |
#btn.setSizePolicy(size_policy) |
#btn.setSizePolicy(size_policy) |
2849 |
2927 |
|
|
2850 |
|
#action = self.toolbar.addAction("redraw", self.recolor) |
|
2851 |
|
#btn = self.toolbar.widgetForAction(action) |
|
2852 |
|
#btn.setAutoRaise(False) |
|
2853 |
|
##btn.setSizePolicy(size_policy) |
|
2854 |
|
|
|
2855 |
|
action = self.toolbar.addAction("hide", self.giracle.hide) |
|
|
2928 |
|
action = self.toolbar.addAction("hide", self.hide) |
2856 |
2929 |
btn = self.toolbar.widgetForAction(action) |
btn = self.toolbar.widgetForAction(action) |
2857 |
2930 |
btn.setAutoRaise(False) |
btn.setAutoRaise(False) |
2858 |
2931 |
#btn.setSizePolicy(size_policy) |
#btn.setSizePolicy(size_policy) |
2859 |
2932 |
|
|
|
2933 |
|
action = self.toolbar.addAction("show", self.show) |
|
2934 |
|
btn = self.toolbar.widgetForAction(action) |
|
2935 |
|
btn.setAutoRaise(False) |
|
2936 |
|
|
|
2937 |
|
action = self.toolbar.addAction("clear", self.clear) |
|
2938 |
|
btn = self.toolbar.widgetForAction(action) |
|
2939 |
|
btn.setAutoRaise(False) |
|
2940 |
|
|
2860 |
2941 |
#self.tool_layout.addWidget(self.toolbar) |
#self.tool_layout.addWidget(self.toolbar) |
2861 |
2942 |
self.addWidget(self.toolbar, row=0, col=0) |
self.addWidget(self.toolbar, row=0, col=0) |
2862 |
2943 |
|
|
|
... |
... |
class HullEstimationWidget(pg.LayoutWidget): |
2901 |
2982 |
self.table = pg.TableWidget(sortable=False) |
self.table = pg.TableWidget(sortable=False) |
2902 |
2983 |
self.splitter.addWidget(self.table) |
self.splitter.addWidget(self.table) |
2903 |
2984 |
|
|
2904 |
|
|
|
2905 |
2985 |
#self.addWidget(self.splitter, row=1, col=0, colspan=5) |
#self.addWidget(self.splitter, row=1, col=0, colspan=5) |
2906 |
2986 |
self.addWidget(self.splitter, row=1, col=0) |
self.addWidget(self.splitter, row=1, col=0) |
|
2987 |
|
|
|
2988 |
|
def hide(self): |
|
2989 |
|
self.giracle.hide() |
|
2990 |
|
self.equation_planes.hide() |
2907 |
2991 |
|
|
|
2992 |
|
def show(self): |
|
2993 |
|
self.giracle.show() |
|
2994 |
|
self.equation_planes.show() |
|
2995 |
|
|
|
2996 |
|
def clear(self): |
|
2997 |
|
self.giracle.clear() |
|
2998 |
|
self.equation_planes.clear() |
|
2999 |
|
|
|
3000 |
|
#оӵ DirectHull понна гинэ |
|
3001 |
|
#č pouze pro DirectHull |
|
3002 |
|
def get_scheme(self): |
|
3003 |
|
scheme = self.param.getValues()['scheme'][0] |
|
3004 |
|
if scheme == 'random': |
|
3005 |
|
ndir = self.param.getValues()['ndir'][0] |
|
3006 |
|
return sball.get_random_directions(ndir, self.ndim) |
|
3007 |
|
elif scheme in self.schemes: |
|
3008 |
|
return self.schemes[scheme].points |
|
3009 |
|
else: |
|
3010 |
|
Scheme = getattr(quadpy.un, scheme) |
|
3011 |
|
self.schemes[scheme] = Scheme(self.ndim) |
|
3012 |
|
return self.schemes[scheme].points |
|
3013 |
|
|
2908 |
3014 |
|
|
2909 |
3015 |
def _set_param(self): |
def _set_param(self): |
2910 |
3016 |
params = list() |
params = list() |
2911 |
|
# params.append({'name': 'method', 'type': 'list', \ |
|
2912 |
|
# 'values': ['Ghull'], 'value': 'Ghull'}) |
|
|
3017 |
|
params.append({'name': 'method', 'type': 'list', 'value': 'DirectHull', \ |
|
3018 |
|
'values': ['SBall', 'BrickHull', 'DirectHull', 'QHull']}) |
|
3019 |
|
params.append({'name': 'space', 'type': 'list', 'tip': "Not used for SBall", \ |
|
3020 |
|
'values': self.sb_item.spaces, 'value': 'G'}) |
|
3021 |
|
|
|
3022 |
|
schemes_list = schemes.un_spheres + ['random'] |
|
3023 |
|
params.append({'name': 'scheme', 'type': 'list', \ |
|
3024 |
|
'values': schemes_list, 'value': schemes_list[0], \ |
|
3025 |
|
'tip': "Used only for DirectHull. Generation can take for a while"}) |
|
3026 |
|
params.append({'name': 'ndir', 'type': 'int', \ |
|
3027 |
|
'limits': (1, float('inf')), 'value': 1000, 'default': 1000, \ |
|
3028 |
|
'title': "number of random directions", \ |
|
3029 |
|
'tip': "Used only for random scheme in DirectHull"}) |
|
3030 |
|
|
2913 |
3031 |
# designs = ['None'] |
# designs = ['None'] |
2914 |
3032 |
# if 'designs' in self.sb_item.kwargs: |
# if 'designs' in self.sb_item.kwargs: |
2915 |
3033 |
# designs.extend(self.sb_item.kwargs['designs'].keys()) |
# designs.extend(self.sb_item.kwargs['designs'].keys()) |
2916 |
3034 |
# params.append({'name': 'design', 'type': 'list', \ |
# params.append({'name': 'design', 'type': 'list', \ |
2917 |
3035 |
# 'values': designs, 'value': 'None'}) |
# 'values': designs, 'value': 'None'}) |
2918 |
3036 |
params.append({'name': 'budget', 'type': 'int', \ |
params.append({'name': 'budget', 'type': 'int', \ |
2919 |
|
'limits': (1, float('inf')), 'value': 1000, 'default': 1000}) |
|
|
3037 |
|
'limits': (1, float('inf')), 'value': 1000, 'default': 1000,\ |
|
3038 |
|
'tip': "Number of simulations for optimal importance sampling"}) |
2920 |
3039 |
|
|
2921 |
3040 |
params.append({'name': 'node (pixel) size', 'type': 'float',\ |
params.append({'name': 'node (pixel) size', 'type': 'float',\ |
2922 |
3041 |
'limits': (0, float('inf')), 'value': 3.5, 'default': self.sb_item.px_size}) |
'limits': (0, float('inf')), 'value': 3.5, 'default': self.sb_item.px_size}) |
2923 |
3042 |
#ё больше калорий богу калорий! |
#ё больше калорий богу калорий! |
2924 |
3043 |
params.append({'name': 'r', 'type': 'color', 'value': (189, 204, 0, 255) }) |
params.append({'name': 'r', 'type': 'color', 'value': (189, 204, 0, 255) }) |
2925 |
3044 |
params.append({'name': 'inside', 'type': 'color', 'value': (133, 172, 102, 255) }) |
params.append({'name': 'inside', 'type': 'color', 'value': (133, 172, 102, 255) }) |
2926 |
|
#params.append({'name': 'convex_hull', 'type': 'color', 'value': (186, 109, 0, 255) }) |
|
|
3045 |
|
params.append({'name': 'convex_hull', 'type': 'color', 'value': (85, 170, 255, 255) }) # (186, 109, 0, 255) |
2927 |
3046 |
params.append({'name': 'fire', 'type': 'color', 'value': (245, 117, 0, 255) }) |
params.append({'name': 'fire', 'type': 'color', 'value': (245, 117, 0, 255) }) |
2928 |
3047 |
params.append({'name': 'outside', 'type': 'color', 'value': 0.6}) |
params.append({'name': 'outside', 'type': 'color', 'value': 0.6}) |
2929 |
|
params.append({'name': 'R', 'type': 'color', 'value': (189, 204, 255, 255) }) |
|
|
3048 |
|
params.append({'name': 'R', 'type': 'color', 'value': (85, 85, 255, 255) }) |
2930 |
3049 |
params.append({'name': 'Update as the box runned', 'type': 'bool', 'value': False }) # 'tip': "This is a checkbox" |
params.append({'name': 'Update as the box runned', 'type': 'bool', 'value': False }) # 'tip': "This is a checkbox" |
2931 |
|
|
|
|
3050 |
|
params.append({'name': 'index', 'title': "replace previous", 'type': 'bool', 'value': True }) |
2932 |
3051 |
|
|
2933 |
3052 |
### Create tree of Parameter objects |
### Create tree of Parameter objects |
2934 |
3053 |
# I don't know why that signals do not work for me |
# I don't know why that signals do not work for me |
|
... |
... |
class HullEstimationWidget(pg.LayoutWidget): |
2937 |
3056 |
#self.param.sigValueChanged.connect(self.param_changed) |
#self.param.sigValueChanged.connect(self.param_changed) |
2938 |
3057 |
#self.param.sigValueChanging.connect(self.param_changing) |
#self.param.sigValueChanging.connect(self.param_changing) |
2939 |
3058 |
self.param = pg.parametertree.Parameter.create(name='params', type='group', children=params) |
self.param = pg.parametertree.Parameter.create(name='params', type='group', children=params) |
|
3059 |
|
|
|
3060 |
|
def get_ghull(self): return khull.Ghull(self.get_hull()) |
|
3061 |
|
#č jistě potřebuji něco, co zpracuje parameter_tree |
|
3062 |
|
#č a vratí platný hull pro Ghull |
|
3063 |
|
def get_hull(self): |
|
3064 |
|
#č bez semplu se neobejde |
|
3065 |
|
nsim = self.sb_item.slider.value() |
|
3066 |
|
sample = self.sb_item.sample_box.f_model[:nsim] |
2940 |
3067 |
|
|
2941 |
|
|
|
2942 |
|
|
|
2943 |
|
|
|
2944 |
|
def recolor(self): |
|
2945 |
|
pass |
|
2946 |
|
# with pg.BusyCursor(): |
|
2947 |
|
# # keep the GUI responsive :) |
|
2948 |
|
# self.sb_item.app.processEvents() |
|
2949 |
|
# #č nejdřív triangulace |
|
2950 |
|
# for tri_bound, plot_item in self.triangulation: |
|
2951 |
|
# plot_item.show() |
|
2952 |
|
# |
|
2953 |
|
# |
|
2954 |
|
# #č teď tečičky |
|
2955 |
|
# for nodes, plot_item, cell_stats in self.simplices: |
|
2956 |
|
# event = cell_stats['event'] |
|
2957 |
|
# if event in self.max_simplices: |
|
2958 |
|
# cell_probability = cell_stats['cell_probability'] |
|
2959 |
|
# cm = self.param.getValues()[event][0] #č očekávám tam kolor mapu |
|
2960 |
|
# blue_intensity = cell_probability / self.max_simplices[event] |
|
2961 |
|
# color = cm.mapToQColor(blue_intensity) |
|
2962 |
|
# else: # outside |
|
2963 |
|
# color = 0.6 |
|
2964 |
|
# #symbolSize = np.sqrt(nodes.w / min(nodes.w)) # not bad |
|
2965 |
|
# size = self.param.getValues()['node (pixel) size'][0] |
|
2966 |
|
# plot_item.setSymbolBrush(color) |
|
2967 |
|
# plot_item.setSymbolSize(size) |
|
2968 |
|
# plot_item.show() |
|
2969 |
|
|
|
2970 |
|
|
|
|
3068 |
|
# ['SBall', 'BrickHull', 'DirectHull', 'QHull'] |
|
3069 |
|
hull_model = self.param.getValues()['method'][0] |
|
3070 |
|
if hull_model == 'SBall': |
|
3071 |
|
return khull.GBall(sample) |
|
3072 |
|
elif hull_model == 'BrickHull': |
|
3073 |
|
space = self.param.getValues()['space'][0] |
|
3074 |
|
return khull.BrickHull(sample, space) |
|
3075 |
|
elif hull_model == 'DirectHull': |
|
3076 |
|
space = self.param.getValues()['space'][0] |
|
3077 |
|
direct_plan = self.get_scheme() |
|
3078 |
|
return khull.DirectHull(sample, direct_plan, space) |
|
3079 |
|
elif hull_model == 'QHull': |
|
3080 |
|
space = self.param.getValues()['space'][0] |
|
3081 |
|
#č tento widget pokažde generuje obálku znovu |
|
3082 |
|
return stm.six.QHull(sample, space, incremental=False) |
|
3083 |
|
else: |
|
3084 |
|
raise ValueError("HullEstimationWidget: co to je za obálku?") |
2971 |
3085 |
|
|
2972 |
3086 |
#č ten hlavní modul se dočkal na překopávání |
#č ten hlavní modul se dočkal na překopávání |
2973 |
3087 |
def on_box_run(self, *args, **kwargs): |
def on_box_run(self, *args, **kwargs): |
2974 |
3088 |
self.giracle.clear() |
self.giracle.clear() |
|
3089 |
|
self.equation_planes.clear() |
2975 |
3090 |
#č je třeba zkontrolovat autorun a restartovat výpočet |
#č je třeba zkontrolovat autorun a restartovat výpočet |
2976 |
3091 |
if self.param.getValues()['Update as the box runned'][0]: |
if self.param.getValues()['Update as the box runned'][0]: |
2977 |
3092 |
self.get_shell_estimation() |
self.get_shell_estimation() |
2978 |
|
#else: |
|
2979 |
|
# self.self_clear() |
|
2980 |
|
|
|
|
3093 |
|
|
|
3094 |
|
def index(self, index): |
|
3095 |
|
if self.param.getValues()['index'][0]: # replace previous |
|
3096 |
|
return index |
|
3097 |
|
else: |
|
3098 |
|
return None |
|
3099 |
|
|
|
3100 |
|
def draw_convex_hull(self, hull): |
|
3101 |
|
try: |
|
3102 |
|
#č zatím uděláme jen pro 2D infinite lajny |
|
3103 |
|
design_points = hull.get_design_points() |
|
3104 |
|
|
|
3105 |
|
if self.param.getValues()['index'][0]: # replace previous |
|
3106 |
|
self.equation_planes.clear() |
|
3107 |
|
|
|
3108 |
|
size = self.param.getValues()['node (pixel) size'][0] |
|
3109 |
|
color = self.param.getValues()['convex_hull'][0] #č tam bude barva |
|
3110 |
|
|
|
3111 |
|
self.giracle.add_serie(design_points, z=31, index=self.index('design points'),\ |
|
3112 |
|
pen=None, symbol='o', symbolPen=pg.mkPen(None), \ |
|
3113 |
|
symbolBrush=color, symbolSize=size, name='design points') |
|
3114 |
|
if self.ndim == 2: |
|
3115 |
|
#č musíme něco zavolat na self.equation_planes |
|
3116 |
|
#č equation_planes má funkci add_line() |
|
3117 |
|
#č add_line(self, space='G', index=None, **plot_kwargs) |
|
3118 |
|
#č která pak plot_kwargs přeposilá funkci addLine() |
|
3119 |
|
#č na central widgetu. |
|
3120 |
|
#č To vše skončí ve pyqtgrafové InfiniteLine třidě. |
|
3121 |
|
#č ta moje třida InfiniteLines sama se stará o shodování prostorů |
|
3122 |
|
#č indexy posilat nebudeme (s nimi je to trošku komplikovanější), |
|
3123 |
|
#č takže ty konvexní obálky budou hromadit |
|
3124 |
|
|
|
3125 |
|
#pos = list() #č navrhové body nakreslíme všechny dohromady |
|
3126 |
|
for equation in hull.equations: |
|
3127 |
|
#č ve 2D bych očekával v rovnici pouze 3 hodnoty (já potřebuji směry) |
|
3128 |
|
x, y, offset = equation |
|
3129 |
|
design_point = [-x*offset, -y*offset] |
|
3130 |
|
#self.sb_item.central_widget.plot(np.array([pos, pos]), symbol='o') |
|
3131 |
|
# if y < 0: #č tak to aspoň kreslí |
|
3132 |
|
# angle = np.rad2deg(np.arcsin(x)) |
|
3133 |
|
# else: |
|
3134 |
|
# angle = np.rad2deg(np.arccos(y)) |
|
3135 |
|
|
|
3136 |
|
if (x*y) < 0: #č tak to aspoň kreslí |
|
3137 |
|
angle = np.rad2deg(np.arccos(np.abs(y))) |
|
3138 |
|
else: |
|
3139 |
|
angle = np.rad2deg(np.arccos(-np.abs(y))) |
|
3140 |
|
self.equation_planes.add_line(space=hull.space,\ |
|
3141 |
|
z=29, pos=design_point, angle=angle, pen=color) |
|
3142 |
|
except BaseException as e: |
|
3143 |
|
msg = "error during estimation " |
|
3144 |
|
error_msg = self.__class__.__name__ + ": " + msg + repr(e) |
|
3145 |
|
print(error_msg) |
|
3146 |
|
self.sb_item.errors.append(e) |
|
3147 |
|
|
|
3148 |
|
def draw_ghull(self, r, R, nodes=None): |
|
3149 |
|
# r-circle |
|
3150 |
|
if r < 0: |
|
3151 |
|
r = 0 |
|
3152 |
|
color = self.param.getValues()['r'][0] #č tam bude barva |
|
3153 |
|
#č krucí, nevím co mám dělat s indexama. |
|
3154 |
|
#č co mám dělat s předchozí kresbou? Nechame |
|
3155 |
|
self.giracle.add_circle(r=r, index=self.index('r'),z=32, pen=color, name='r') |
|
3156 |
|
|
|
3157 |
|
# R-circle. #č Při kreslení nahradíme předchozí. |
|
3158 |
|
color = self.param.getValues()['R'][0] #č tam bude barva |
|
3159 |
|
self.giracle.add_circle(R, z=32, index=self.index('R'), pen=color, name='R') |
|
3160 |
|
|
|
3161 |
|
#č draw tečičky |
|
3162 |
|
if nodes is not None: |
|
3163 |
|
size = self.param.getValues()['node (pixel) size'][0] |
|
3164 |
|
plot_params = {'pen':None, 'symbol':'o', 'symbolSize':size, \ |
|
3165 |
|
'symbolPen':pg.mkPen(None)} |
|
3166 |
|
#brush = pg.mkBrush(color) |
|
3167 |
|
mask = nodes.is_outside |
|
3168 |
|
|
|
3169 |
|
index = 'inside' |
|
3170 |
|
color = self.param.getValues()[index][0] #č tam bude barva |
|
3171 |
|
self.giracle.add_serie(nodes[~mask], z=30, index=self.index(index),\ |
|
3172 |
|
symbolBrush=color, name=index, **plot_params) |
|
3173 |
|
|
|
3174 |
|
index = 'outside' |
|
3175 |
|
color = self.param.getValues()[index][0] #č tam bude barva |
|
3176 |
|
self.giracle.add_serie(nodes[mask], z=30, index=self.index(index),\ |
|
3177 |
|
symbolBrush=color, name=index, **plot_params) |
2981 |
3178 |
|
|
2982 |
|
# def self_clear(self): |
|
2983 |
|
# # odebereme prvky-propísky z hlavního plotu |
|
2984 |
|
# for tri_bound, plot_item in self.triangulation: |
|
2985 |
|
# self.sb_item.central_widget.removeItem(plot_item) |
|
2986 |
|
# for nodes, plot_item, cell_stats in self.simplices: |
|
2987 |
|
# self.sb_item.central_widget.removeItem(plot_item) |
|
2988 |
|
# |
|
2989 |
|
# self.redraw() |
|
2990 |
3179 |
|
|
2991 |
3180 |
def get_shell_estimation(self): |
def get_shell_estimation(self): |
2992 |
|
nsim = self.sb_item.slider.value() |
|
2993 |
|
sample = self.sb_item.sample_box.f_model[:nsim] |
|
|
3181 |
|
ghull = self.get_ghull() |
|
3182 |
|
self.draw_convex_hull(ghull.hull) |
2994 |
3183 |
|
|
2995 |
3184 |
try: |
try: |
2996 |
|
ghull = stm.six.Ghull(sample, incremental=False, design=None) |
|
2997 |
3185 |
shell_estimation, global_stats = ghull.get_shell_estimation() |
shell_estimation, global_stats = ghull.get_shell_estimation() |
2998 |
3186 |
|
|
2999 |
3187 |
# if hasattr(self.sb_item.sample_box, 'estimations'): |
# if hasattr(self.sb_item.sample_box, 'estimations'): |
|
... |
... |
class HullEstimationWidget(pg.LayoutWidget): |
3001 |
3189 |
# self.sb_item.estimation_added.emit() |
# self.sb_item.estimation_added.emit() |
3002 |
3190 |
self.table.setData({**global_stats, 'shell_estimation':shell_estimation}) |
self.table.setData({**global_stats, 'shell_estimation':shell_estimation}) |
3003 |
3191 |
|
|
3004 |
|
index = 'r' |
|
3005 |
|
r = global_stats[index] |
|
3006 |
|
if r < 0: |
|
3007 |
|
r = 0 |
|
3008 |
|
color = self.param.getValues()[index][0] #č bude tam barva? |
|
3009 |
|
plot_item = self.giracle.add_circle(r=r,\ |
|
3010 |
|
index=index, pen=color, name=index) |
|
3011 |
|
plot_item.setZValue(32) |
|
3012 |
|
|
|
3013 |
|
index = 'R' |
|
3014 |
|
color = self.param.getValues()[index][0] #č bude tam barva? |
|
3015 |
|
plot_item = self.giracle.add_circle(global_stats[index],\ |
|
3016 |
|
index=index, pen=color, name=index) |
|
3017 |
|
plot_item.setZValue(32) |
|
|
3192 |
|
self.draw_ghull(global_stats['r'], global_stats['R']) |
3018 |
3193 |
|
|
3019 |
3194 |
except BaseException as e: |
except BaseException as e: |
3020 |
3195 |
msg = "error during estimation " |
msg = "error during estimation " |
3021 |
3196 |
error_msg = self.__class__.__name__ + ": " + msg + repr(e) |
error_msg = self.__class__.__name__ + ": " + msg + repr(e) |
3022 |
3197 |
print(error_msg) |
print(error_msg) |
|
3198 |
|
self.sb_item.errors.append(e) |
3023 |
3199 |
|
|
3024 |
3200 |
|
|
3025 |
3201 |
def fire(self): |
def fire(self): |
3026 |
|
nsim = self.sb_item.slider.value() |
|
3027 |
|
sample = self.sb_item.sample_box.f_model[:nsim] |
|
|
3202 |
|
#č zatím máme issue s náhodným planem |
|
3203 |
|
#č kdyby něco - vyřeším přes ukladání |
|
3204 |
|
#č náhodných plánů v parameter_tree |
|
3205 |
|
hull = self.get_hull() |
|
3206 |
|
self.draw_convex_hull(hull) |
3028 |
3207 |
#☺ Krucinal, kdo ten OrderedDict vymyslel? |
#☺ Krucinal, kdo ten OrderedDict vymyslel? |
3029 |
3208 |
params = self.param.getValues() |
params = self.param.getValues() |
3030 |
3209 |
ns = params['budget'][0] |
ns = params['budget'][0] |
3031 |
3210 |
try: |
try: |
3032 |
|
ghull = stm.six.Ghull(sample, incremental=False, design=None) |
|
3033 |
|
fire_G = ghull.fire(ns) |
|
3034 |
|
fire = sample.new_sample(fire_G, space='G', extend=True) |
|
3035 |
|
|
|
3036 |
|
color = self.param.getValues()['fire'][0] #č bude tam barva? |
|
3037 |
|
|
|
3038 |
|
|
|
|
3211 |
|
fire = hull.fire(ns) |
3039 |
3212 |
# draw tečičky |
# draw tečičky |
3040 |
|
|
|
|
3213 |
|
color = self.param.getValues()['fire'][0] #č tam bude barva |
3041 |
3214 |
size = self.param.getValues()['node (pixel) size'][0] |
size = self.param.getValues()['node (pixel) size'][0] |
3042 |
3215 |
#brush = pg.mkBrush(color) |
#brush = pg.mkBrush(color) |
3043 |
|
plot_item = self.giracle.add_serie(fire, index='fire',\ |
|
|
3216 |
|
self.giracle.add_serie(fire, z=30, index=self.index('fire'),\ |
3044 |
3217 |
pen=None, symbol='o', symbolPen=pg.mkPen(None), \ |
pen=None, symbol='o', symbolPen=pg.mkPen(None), \ |
3045 |
3218 |
symbolBrush=color, symbolSize=size, name='fire') |
symbolBrush=color, symbolSize=size, name='fire') |
3046 |
|
plot_item.setZValue(30) |
|
3047 |
3219 |
|
|
3048 |
3220 |
except BaseException as e: |
except BaseException as e: |
3049 |
3221 |
msg = "error " |
msg = "error " |
3050 |
3222 |
error_msg = self.__class__.__name__ + ": " + msg + repr(e) |
error_msg = self.__class__.__name__ + ": " + msg + repr(e) |
3051 |
3223 |
print(error_msg) |
print(error_msg) |
3052 |
|
|
|
|
3224 |
|
self.sb_item.errors.append(e) |
3053 |
3225 |
|
|
3054 |
3226 |
def integrate(self): |
def integrate(self): |
3055 |
|
nsim = self.sb_item.slider.value() |
|
3056 |
|
sample = self.sb_item.sample_box.f_model[:nsim] |
|
|
3227 |
|
ghull = self.get_ghull() |
|
3228 |
|
self.draw_convex_hull(ghull.hull) |
|
3229 |
|
|
3057 |
3230 |
#☺ Krucinal, kdo ten OrderedDict vymyslel? |
#☺ Krucinal, kdo ten OrderedDict vymyslel? |
3058 |
3231 |
params = self.param.getValues() |
params = self.param.getValues() |
3059 |
3232 |
budget = params['budget'][0] |
budget = params['budget'][0] |
|
... |
... |
class HullEstimationWidget(pg.LayoutWidget): |
3065 |
3238 |
# design = self.sb_item.kwargs['designs'][design] |
# design = self.sb_item.kwargs['designs'][design] |
3066 |
3239 |
|
|
3067 |
3240 |
try: |
try: |
3068 |
|
ghull = stm.six.Ghull(sample, incremental=False, design=None) |
|
3069 |
3241 |
data = ghull.integrate(budget) |
data = ghull.integrate(budget) |
3070 |
3242 |
nodes, ghull_estimation, convex_hull_estimation, global_stats = data |
nodes, ghull_estimation, convex_hull_estimation, global_stats = data |
3071 |
3243 |
|
|
|
... |
... |
class HullEstimationWidget(pg.LayoutWidget): |
3075 |
3247 |
self.table.setData({**global_stats, "ghull_estimation":ghull_estimation,\ |
self.table.setData({**global_stats, "ghull_estimation":ghull_estimation,\ |
3076 |
3248 |
"convex_hull_estimation": convex_hull_estimation}) |
"convex_hull_estimation": convex_hull_estimation}) |
3077 |
3249 |
|
|
3078 |
|
# draw tečičky |
|
3079 |
|
size = self.param.getValues()['node (pixel) size'][0] |
|
3080 |
|
plot_params = {'pen':None, 'symbol':'o', 'symbolSize':size, \ |
|
3081 |
|
'symbolPen':pg.mkPen(None)} |
|
3082 |
|
#brush = pg.mkBrush(color) |
|
3083 |
|
mask = nodes.is_outside |
|
3084 |
|
|
|
3085 |
|
index = 'inside' |
|
3086 |
|
color = self.param.getValues()[index][0] #č bude tam barva? |
|
3087 |
|
plot_item = self.giracle.add_serie(nodes[mask], index=index,\ |
|
3088 |
|
symbolBrush=color, name=index, **plot_params) |
|
3089 |
|
plot_item.setZValue(30) |
|
3090 |
|
|
|
3091 |
|
index = 'outside' |
|
3092 |
|
color = self.param.getValues()[index][0] #č bude tam barva? |
|
3093 |
|
plot_item = self.giracle.add_serie(nodes[~mask], index=index,\ |
|
3094 |
|
symbolBrush=color, name=index, **plot_params) |
|
3095 |
|
plot_item.setZValue(30) |
|
3096 |
|
|
|
3097 |
|
index = 'r' |
|
3098 |
|
r = global_stats[index] |
|
3099 |
|
if r < 0: |
|
3100 |
|
r = 0 |
|
3101 |
|
color = self.param.getValues()[index][0] #č bude tam barva? |
|
3102 |
|
plot_item = self.giracle.add_circle(r=r,\ |
|
3103 |
|
index=index, pen=color, name=index) |
|
3104 |
|
plot_item.setZValue(32) |
|
3105 |
|
|
|
3106 |
|
index = 'R' |
|
3107 |
|
color = self.param.getValues()[index][0] #č bude tam barva? |
|
3108 |
|
plot_item = self.giracle.add_circle(global_stats[index],\ |
|
3109 |
|
index=index, pen=color, name=index) |
|
3110 |
|
plot_item.setZValue(32) |
|
|
3250 |
|
self.draw_ghull(global_stats['r'], global_stats['R'], nodes) |
3111 |
3251 |
|
|
3112 |
3252 |
|
|
3113 |
3253 |
except BaseException as e: |
except BaseException as e: |
3114 |
3254 |
msg = "error during estimation " |
msg = "error during estimation " |
3115 |
3255 |
error_msg = self.__class__.__name__ + ": " + msg + repr(e) |
error_msg = self.__class__.__name__ + ": " + msg + repr(e) |
3116 |
3256 |
print(error_msg) |
print(error_msg) |
3117 |
|
|
|
|
3257 |
|
self.sb_item.errors.append(e) |
3118 |
3258 |
|
|
3119 |
3259 |
|
|
3120 |
3260 |
|
|
File simplex.py changed (mode: 100644) (index 9c0f513..c804dbc) |
... |
... |
from scipy import spatial |
9 |
9 |
from scipy import stats |
from scipy import stats |
10 |
10 |
from .candybox import CandyBox |
from .candybox import CandyBox |
11 |
11 |
|
|
12 |
|
#č mým úkolem při návrhu této třidy je pořádně všecko zkomplikovat. |
|
13 |
|
#č Dostávám za to peníze. |
|
14 |
|
#č Takže. Udělám 3 druhů estimátorů |
|
15 |
|
# convex_hull_estimation -2: inside, -1: outside |
|
16 |
|
# shell_estimation -22: inner, -3: shell, -11: outer |
|
17 |
|
# ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer |
|
18 |
|
class Ghull: |
|
19 |
|
def __init__(self, sample, incremental=True, design=None): |
|
|
12 |
|
#č jako sůl potřebujem statefull třidu, |
|
13 |
|
#č která by schovávala vnitřní implementaciju, |
|
14 |
|
#č která by nabízela jednotné rozhraní pro ostatní kód WellMet |
|
15 |
|
#č a byla by přiměřeně kompatibilní s ConvexHull ze scipy.spatial |
|
16 |
|
#č Ze scipy bych viděl podporu atributů .points, .npoints, .equations |
|
17 |
|
#č Určitě bych chtěl funkce .is_outside(sample) |
|
18 |
|
#č Atribute .space by ukazoval, v jakém prostoru konvexní obál je vytvořen |
|
19 |
|
|
|
20 |
|
|
|
21 |
|
class QHull: |
|
22 |
|
def __init__(self, sample, space='G', incremental=True): |
20 |
23 |
self.sample = sample |
self.sample = sample |
21 |
|
self.design = design |
|
22 |
|
self.convex_hull = spatial.ConvexHull(sample.G, incremental=incremental) |
|
23 |
|
self.shell = sball.Shell(sample.nvar) |
|
24 |
|
|
|
25 |
|
self.calculate_d() |
|
26 |
|
|
|
27 |
|
def update(self): |
|
28 |
|
self.convex_hull.add_points(self.sample[self.convex_hull.npoints:].G) |
|
29 |
|
self.calculate_d() |
|
30 |
|
|
|
31 |
|
def get_R(self): |
|
32 |
|
sum_squared = np.sum(np.square(self.sample.G), axis=1) |
|
33 |
|
#index = np.argmax(sum_squared) |
|
34 |
|
return np.sqrt(np.nanmax(sum_squared)) |
|
|
24 |
|
self.incremental = incremental |
|
25 |
|
self.space = space |
|
26 |
|
|
|
27 |
|
def regen(self): |
|
28 |
|
points = getattr(self.sample, self.space) |
|
29 |
|
self.convex_hull = spatial.ConvexHull(points, incremental=self.incremental) |
|
30 |
|
|
|
31 |
|
def __getattr__(self, attr): |
|
32 |
|
#č branime se rekurzii |
|
33 |
|
# defend against recursion |
|
34 |
|
#оӵ рекурсилы пезьдэт! |
|
35 |
|
if attr == 'convex_hull': |
|
36 |
|
#č zkusme rychle konvexní obálky sestavit |
|
37 |
|
#č a ihned ji vrátit |
|
38 |
|
self.regen() |
|
39 |
|
return self.convex_hull |
|
40 |
|
|
|
41 |
|
elif attr == 'A': |
|
42 |
|
return self.convex_hull.equations[:,:-1] |
|
43 |
|
elif attr == 'b': |
|
44 |
|
return self.convex_hull.equations[:,-1] |
|
45 |
|
|
|
46 |
|
#ё По всем вопросам обращайтесь |
|
47 |
|
#ё на нашу горячую линию |
|
48 |
|
else: |
|
49 |
|
self._update #č dycky čerstý chleba! |
|
50 |
|
return getattr(self.convex_hull, attr) |
|
51 |
|
|
|
52 |
|
|
|
53 |
|
#č nejsem jist, jestli ten update vůbec dělat. |
|
54 |
|
#č lze navrhnout třidu aj tak, že sama bude hlídat změny. |
|
55 |
|
#č Jenomže, co kdybychom ten automatický update nechtěli? |
|
56 |
|
def _update(self): |
|
57 |
|
if self.convex_hull.npoints < self.sample.nsim: |
|
58 |
|
if self.incremental: |
|
59 |
|
points = getattr(self.sample, self.space) |
|
60 |
|
self.convex_hull.add_points(points[self.convex_hull.npoints:]) |
|
61 |
|
else: |
|
62 |
|
self.regen() |
35 |
63 |
|
|
36 |
|
def calculate_d(self): |
|
37 |
|
"calculates distances to planes. Can be negative" |
|
38 |
|
A = self.convex_hull.equations[:,:-1] |
|
39 |
|
b = self.convex_hull.equations[:,-1] |
|
40 |
|
sum_squared = np.sum(np.square(A), axis=1) |
|
41 |
|
self.d = b / np.sqrt(sum_squared) |
|
42 |
64 |
|
|
43 |
|
def get_r(self): |
|
44 |
|
return -np.nanmax(self.d) |
|
|
65 |
|
def is_inside(self, nodes): |
|
66 |
|
self._update() |
|
67 |
|
x = getattr(nodes, self.space) |
45 |
68 |
|
|
46 |
|
def fire(self, ns): |
|
|
69 |
|
#E [normal, offset] forming the hyperplane equation of the facet (see Qhull documentation for more) |
47 |
70 |
A = self.convex_hull.equations[:,:-1] |
A = self.convex_hull.equations[:,:-1] |
48 |
71 |
b = self.convex_hull.equations[:,-1] |
b = self.convex_hull.equations[:,-1] |
49 |
72 |
|
|
50 |
|
to_fire = np.argmax(self.d) |
|
51 |
|
a = A[to_fire] |
|
52 |
|
fire_from = stats.norm.cdf(self.get_r()) |
|
53 |
|
t = np.linspace(fire_from, 1, ns) |
|
54 |
|
t = stats.norm.ppf(t) |
|
55 |
|
fire_G = t.reshape(-1,1) @ a.reshape(1,-1) |
|
|
73 |
|
# N=ns, E - number of hyperplane equations |
|
74 |
|
ExN = A @ x.T + np.atleast_2d(b).T |
|
75 |
|
mask = np.all(ExN < 0, axis=0) |
|
76 |
|
return mask |
|
77 |
|
|
|
78 |
|
def is_outside(hull, nodes): |
|
79 |
|
return ~hull.is_inside(nodes) |
56 |
80 |
|
|
57 |
|
return fire_G |
|
|
81 |
|
def get_design_points(hull): |
|
82 |
|
hull._update() |
|
83 |
|
sample_model = -hull.A * hull.b.reshape(-1,1) |
|
84 |
|
return hull.sample.f_model.new_sample(sample_model, space=hull.space) |
58 |
85 |
|
|
59 |
|
def get_shell_estimation(self): |
|
60 |
|
shell = self.shell |
|
61 |
|
r = self.get_r() |
|
62 |
|
R = self.get_R() |
|
63 |
86 |
|
|
64 |
|
if r<0: |
|
65 |
|
shell.set_bounds(0, R) |
|
|
87 |
|
# shortcut for Ghull |
|
88 |
|
# valid only if space==G |
|
89 |
|
def get_r(hull): |
|
90 |
|
if hull.space=='G': |
|
91 |
|
hull._update() |
|
92 |
|
b = hull.convex_hull.equations[:,-1] |
|
93 |
|
return -np.nanmax(b) |
66 |
94 |
else: |
else: |
67 |
|
shell.set_bounds(r, R) |
|
68 |
|
|
|
69 |
|
# shell_estimation -22: inner, -3: shell, -11: outer |
|
70 |
|
shell_estimation = {-22:shell.ps, -3: shell.p_shell, -11: shell.pf} |
|
71 |
|
global_stats = {"r":r, "R":R, "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf} |
|
72 |
|
return shell_estimation, global_stats |
|
73 |
|
|
|
74 |
|
def integrate(self, nis): |
|
75 |
|
#č no to teda disajn třidy je doopravdy hroznej |
|
76 |
|
# .get_shell_estimation() will calculate radiuses and will update shell |
|
77 |
|
shell_estimation, global_stats = self.get_shell_estimation() |
|
78 |
|
nodes_G = self.shell.rvs(nis) |
|
79 |
|
mask = is_outside(self.convex_hull, nodes_G) |
|
80 |
|
#č nevím proč, ale kdysi mě to vyšlo rychlejc jak obyčejný součet |
|
81 |
|
nf = len(mask[mask]) |
|
82 |
|
ns = nis - nf |
|
83 |
|
# shell_estimation -22: inner, -3: shell, -11: outer |
|
84 |
|
p_shell = shell_estimation.pop(-3) |
|
85 |
|
shell_pf = nf/nis * p_shell |
|
86 |
|
shell_ps = ns/nis * p_shell |
|
87 |
|
|
|
88 |
|
# ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer |
|
89 |
|
ghull_estimation = shell_estimation; del(shell_estimation) |
|
90 |
|
ghull_estimation[-21] = shell_ps |
|
91 |
|
ghull_estimation[-12] = shell_pf |
|
92 |
|
global_stats["shell inside"] = shell_ps |
|
93 |
|
global_stats["shell outside"] = shell_pf |
|
94 |
|
|
|
95 |
|
# convex_hull_estimation -2: inside, -1: outside |
|
96 |
|
inside = shell_ps + self.shell.ps |
|
97 |
|
outside = shell_pf + self.shell.pf |
|
98 |
|
convex_hull_estimation = {-2: inside, -1: outside} |
|
99 |
|
global_stats["inside"] = inside |
|
100 |
|
global_stats["outside"] = outside |
|
101 |
|
|
|
102 |
|
nodes = self.sample.f_model.new_sample(nodes_G, space='G') |
|
103 |
|
# -2 = 'inside' -1 = 'outside' |
|
104 |
|
nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask) |
|
105 |
|
|
|
106 |
|
return nodes, ghull_estimation, convex_hull_estimation, global_stats |
|
|
95 |
|
return 0 |
|
96 |
|
|
|
97 |
|
def fire(hull, ns): |
|
98 |
|
if hull.space == 'G': |
|
99 |
|
A = hull.equations[:,:-1] |
|
100 |
|
b = hull.equations[:,-1] |
|
101 |
|
|
|
102 |
|
to_fire = np.argmax(b) |
|
103 |
|
a = A[to_fire] |
|
104 |
|
fire_from = stats.norm.cdf(hull.get_r()) |
|
105 |
|
t = np.linspace(fire_from, 1, ns) |
|
106 |
|
t = stats.norm.ppf(t) |
|
107 |
|
fire_G = t.reshape(-1,1) @ a.reshape(1,-1) |
|
108 |
|
|
|
109 |
|
return hull.sample.f_model.new_sample(fire_G, space='G') |
|
110 |
|
|
|
111 |
|
##č mým úkolem při návrhu této třidy je pořádně všecko zkomplikovat. |
|
112 |
|
##č Dostávám za to peníze. |
|
113 |
|
##č Takže. Udělám 3 druhů estimátorů |
|
114 |
|
## convex_hull_estimation -2: inside, -1: outside |
|
115 |
|
## shell_estimation -22: inner, -3: shell, -11: outer |
|
116 |
|
## ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer |
|
117 |
|
#class Ghull: |
|
118 |
|
# def __init__(self, hull, design=None): |
|
119 |
|
# self.hull = hull |
|
120 |
|
# self.design = design |
|
121 |
|
# self.shell = sball.Shell(hull.sample.nvar) |
|
122 |
|
# self.sample = hull.sample |
|
123 |
|
# |
|
124 |
|
# self.fire = hull.fire |
|
125 |
|
# |
|
126 |
|
# def get_R(self): |
|
127 |
|
# sum_squared = np.sum(np.square(self.sample.G), axis=1) |
|
128 |
|
# #index = np.argmax(sum_squared) |
|
129 |
|
# return np.sqrt(np.nanmax(sum_squared)) |
|
130 |
|
# |
|
131 |
|
# def get_r(self): |
|
132 |
|
# "calculates minimal signed distance to planes. Can therefore be negative" |
|
133 |
|
# #return -np.nanmax(self.hull.b) |
|
134 |
|
# return self.hull.get_r() |
|
135 |
|
# |
|
136 |
|
# def get_shell_estimation(self): |
|
137 |
|
# shell = self.shell |
|
138 |
|
# r = self.get_r() |
|
139 |
|
# R = self.get_R() |
|
140 |
|
# |
|
141 |
|
# if r<0: |
|
142 |
|
# shell.set_bounds(0, R) |
|
143 |
|
# else: |
|
144 |
|
# shell.set_bounds(r, R) |
|
145 |
|
# |
|
146 |
|
# # shell_estimation -22: inner, -3: shell, -11: outer |
|
147 |
|
# shell_estimation = {-22:shell.ps, -3: shell.p_shell, -11: shell.pf} |
|
148 |
|
# global_stats = {"r":r, "R":R, "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf} |
|
149 |
|
# return shell_estimation, global_stats |
|
150 |
|
# |
|
151 |
|
# def integrate(self, nis): |
|
152 |
|
# #č no to teda disajn třidy je doopravdy hroznej |
|
153 |
|
# # .get_shell_estimation() will calculate radiuses and will update shell |
|
154 |
|
# shell_estimation, global_stats = self.get_shell_estimation() |
|
155 |
|
# nodes_G = self.shell.rvs(nis) |
|
156 |
|
# nodes = self.sample.f_model.new_sample(nodes_G, space='G') |
|
157 |
|
# mask = self.hull.is_outside(nodes) |
|
158 |
|
# #č nevím proč, ale kdysi mě to vyšlo rychlejc jak obyčejný součet |
|
159 |
|
# nf = len(mask[mask]) |
|
160 |
|
# ns = nis - nf |
|
161 |
|
# # shell_estimation -22: inner, -3: shell, -11: outer |
|
162 |
|
# p_shell = shell_estimation.pop(-3) |
|
163 |
|
# shell_pf = nf/nis * p_shell |
|
164 |
|
# shell_ps = ns/nis * p_shell |
|
165 |
|
# |
|
166 |
|
# # ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer |
|
167 |
|
# ghull_estimation = shell_estimation; del(shell_estimation) |
|
168 |
|
# ghull_estimation[-21] = shell_ps |
|
169 |
|
# ghull_estimation[-12] = shell_pf |
|
170 |
|
# global_stats["shell inside"] = shell_ps |
|
171 |
|
# global_stats["shell outside"] = shell_pf |
|
172 |
|
# |
|
173 |
|
# # convex_hull_estimation -2: inside, -1: outside |
|
174 |
|
# inside = shell_ps + self.shell.ps |
|
175 |
|
# outside = shell_pf + self.shell.pf |
|
176 |
|
# convex_hull_estimation = {-2: inside, -1: outside} |
|
177 |
|
# global_stats["inside"] = inside |
|
178 |
|
# global_stats["outside"] = outside |
|
179 |
|
# |
|
180 |
|
# #nodes = self.sample.f_model.new_sample(nodes_G, space='G') |
|
181 |
|
# # -2 = 'inside' -1 = 'outside' |
|
182 |
|
# nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask) |
|
183 |
|
# |
|
184 |
|
# return nodes, ghull_estimation, convex_hull_estimation, global_stats |
|
185 |
|
# |
|
186 |
|
# |
|
187 |
|
##č mým úkolem při návrhu této třidy je pořádně všecko zkomplikovat. |
|
188 |
|
##č Dostávám za to peníze. |
|
189 |
|
##č Takže. Udělám 3 druhů estimátorů |
|
190 |
|
## convex_hull_estimation -2: inside, -1: outside |
|
191 |
|
## shell_estimation -22: inner, -3: shell, -11: outer |
|
192 |
|
## ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer |
|
193 |
|
#class Ghull: |
|
194 |
|
# def __init__(self, sample, incremental=True, design=None): |
|
195 |
|
# self.sample = sample |
|
196 |
|
# self.design = design |
|
197 |
|
# self.convex_hull = spatial.ConvexHull(sample.G, incremental=incremental) |
|
198 |
|
# self.shell = sball.Shell(sample.nvar) |
|
199 |
|
# |
|
200 |
|
# self.calculate_d() |
|
201 |
|
# |
|
202 |
|
# def update(self): |
|
203 |
|
# self.convex_hull.add_points(self.sample[self.convex_hull.npoints:].G) |
|
204 |
|
# self.calculate_d() |
|
205 |
|
# |
|
206 |
|
# def get_R(self): |
|
207 |
|
# sum_squared = np.sum(np.square(self.sample.G), axis=1) |
|
208 |
|
# #index = np.argmax(sum_squared) |
|
209 |
|
# return np.sqrt(np.nanmax(sum_squared)) |
|
210 |
|
# |
|
211 |
|
# def calculate_d(self): |
|
212 |
|
# "calculates distances to planes. Can be negative" |
|
213 |
|
# A = self.convex_hull.equations[:,:-1] |
|
214 |
|
# b = self.convex_hull.equations[:,-1] |
|
215 |
|
# sum_squared = np.sum(np.square(A), axis=1) |
|
216 |
|
# self.d = b / np.sqrt(sum_squared) |
|
217 |
|
# |
|
218 |
|
# def get_r(self): |
|
219 |
|
# return -np.nanmax(self.d) |
|
220 |
|
# |
|
221 |
|
# def fire(self, ns): |
|
222 |
|
# A = self.convex_hull.equations[:,:-1] |
|
223 |
|
# b = self.convex_hull.equations[:,-1] |
|
224 |
|
# |
|
225 |
|
# to_fire = np.argmax(self.d) |
|
226 |
|
# a = A[to_fire] |
|
227 |
|
# fire_from = stats.norm.cdf(self.get_r()) |
|
228 |
|
# t = np.linspace(fire_from, 1, ns) |
|
229 |
|
# t = stats.norm.ppf(t) |
|
230 |
|
# fire_G = t.reshape(-1,1) @ a.reshape(1,-1) |
|
231 |
|
# |
|
232 |
|
# return fire_G |
|
233 |
|
# |
|
234 |
|
# def get_shell_estimation(self): |
|
235 |
|
# shell = self.shell |
|
236 |
|
# r = self.get_r() |
|
237 |
|
# R = self.get_R() |
|
238 |
|
# |
|
239 |
|
# if r<0: |
|
240 |
|
# shell.set_bounds(0, R) |
|
241 |
|
# else: |
|
242 |
|
# shell.set_bounds(r, R) |
|
243 |
|
# |
|
244 |
|
# # shell_estimation -22: inner, -3: shell, -11: outer |
|
245 |
|
# shell_estimation = {-22:shell.ps, -3: shell.p_shell, -11: shell.pf} |
|
246 |
|
# global_stats = {"r":r, "R":R, "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf} |
|
247 |
|
# return shell_estimation, global_stats |
|
248 |
|
# |
|
249 |
|
# def integrate(self, nis): |
|
250 |
|
# #č no to teda disajn třidy je doopravdy hroznej |
|
251 |
|
# # .get_shell_estimation() will calculate radiuses and will update shell |
|
252 |
|
# shell_estimation, global_stats = self.get_shell_estimation() |
|
253 |
|
# nodes_G = self.shell.rvs(nis) |
|
254 |
|
# mask = is_outside(self.convex_hull, nodes_G) |
|
255 |
|
# #č nevím proč, ale kdysi mě to vyšlo rychlejc jak obyčejný součet |
|
256 |
|
# nf = len(mask[mask]) |
|
257 |
|
# ns = nis - nf |
|
258 |
|
# # shell_estimation -22: inner, -3: shell, -11: outer |
|
259 |
|
# p_shell = shell_estimation.pop(-3) |
|
260 |
|
# shell_pf = nf/nis * p_shell |
|
261 |
|
# shell_ps = ns/nis * p_shell |
|
262 |
|
# |
|
263 |
|
# # ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer |
|
264 |
|
# ghull_estimation = shell_estimation; del(shell_estimation) |
|
265 |
|
# ghull_estimation[-21] = shell_ps |
|
266 |
|
# ghull_estimation[-12] = shell_pf |
|
267 |
|
# global_stats["shell inside"] = shell_ps |
|
268 |
|
# global_stats["shell outside"] = shell_pf |
|
269 |
|
# |
|
270 |
|
# # convex_hull_estimation -2: inside, -1: outside |
|
271 |
|
# inside = shell_ps + self.shell.ps |
|
272 |
|
# outside = shell_pf + self.shell.pf |
|
273 |
|
# convex_hull_estimation = {-2: inside, -1: outside} |
|
274 |
|
# global_stats["inside"] = inside |
|
275 |
|
# global_stats["outside"] = outside |
|
276 |
|
# |
|
277 |
|
# nodes = self.sample.f_model.new_sample(nodes_G, space='G') |
|
278 |
|
# # -2 = 'inside' -1 = 'outside' |
|
279 |
|
# nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask) |
|
280 |
|
# |
|
281 |
|
# return nodes, ghull_estimation, convex_hull_estimation, global_stats |
107 |
282 |
|
|
108 |
283 |
|
|
109 |
284 |
|
|