File simplex.py changed (mode: 100644) (index c804dbc..b41bdcc) |
... |
... |
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 |
|
#č 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): |
|
23 |
|
self.sample = sample |
|
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() |
|
63 |
|
|
|
64 |
|
|
|
65 |
|
def is_inside(self, nodes): |
|
66 |
|
self._update() |
|
67 |
|
x = getattr(nodes, self.space) |
|
68 |
|
|
|
69 |
|
#E [normal, offset] forming the hyperplane equation of the facet (see Qhull documentation for more) |
|
70 |
|
A = self.convex_hull.equations[:,:-1] |
|
71 |
|
b = self.convex_hull.equations[:,-1] |
|
72 |
|
|
|
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) |
|
80 |
|
|
|
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) |
|
85 |
|
|
|
86 |
|
|
|
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) |
|
94 |
|
else: |
|
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 |
|
282 |
|
|
|
283 |
|
|
|
284 |
|
|
|
285 |
|
|
|
286 |
12 |
#č napadlo mě zababáchnuť třidu, která by se sama starala o všem co se tyče |
#č napadlo mě zababáchnuť třidu, která by se sama starala o všem co se tyče |
287 |
13 |
#č vnější domény. Nešlo mě totíž to udělat jednou funkcí, bylo by velmi |
#č vnější domény. Nešlo mě totíž to udělat jednou funkcí, bylo by velmi |
288 |
14 |
#č špatné z hlediska zodpovednosti kódu. Tak to všecko zabalíme to třidy |
#č špatné z hlediska zodpovednosti kódu. Tak to všecko zabalíme to třidy |