File wellmet/voronoi.py changed (mode: 100644) (index 331b9e2..04923de) |
... |
... |
class ConvexSolver: |
201 |
201 |
|
|
202 |
202 |
|
|
203 |
203 |
|
|
|
204 |
|
#č vlastní slovník se dycky hodí |
|
205 |
|
class WardrobeDict(dict): |
|
206 |
|
def __init__(self): |
|
207 |
|
super().__init__() |
|
208 |
|
self.index = 0 |
|
209 |
|
|
|
210 |
|
def add(self, item): |
|
211 |
|
self.index += 1 |
|
212 |
|
#оӵ мар ке вал отын - туж жаль |
|
213 |
|
self[self.index] = item |
|
214 |
|
return self.index |
204 |
215 |
|
|
205 |
216 |
|
|
206 |
217 |
#оӵ кык точкаен Вороной |
#оӵ кык точкаен Вороной |
|
... |
... |
class ContactVoronoi: |
240 |
251 |
č a mezi "c" a "b", tak přídaní "c" nemohlo kontakt mezi "a" a "b" ovlivnit |
č a mezi "c" a "b", tak přídaní "c" nemohlo kontakt mezi "a" a "b" ovlivnit |
241 |
252 |
""" |
""" |
242 |
253 |
def __init__(self, sample_box, hull, model_space='G', \ |
def __init__(self, sample_box, hull, model_space='G', \ |
243 |
|
p_norm=2, p_base=0.5, ns=1000, auto_update=True, \ |
|
244 |
|
sample_callback=None): |
|
|
254 |
|
p_norm=2, ns=1000, p_base=0.5, auto_update=True, \ |
|
255 |
|
on_add_mixed=None, on_update_mixed=None, on_delete_mixed=None): |
245 |
256 |
self.sample_box = sample_box |
self.sample_box = sample_box |
246 |
257 |
self.hull = hull |
self.hull = hull |
247 |
258 |
self.model_space = model_space |
self.model_space = model_space |
248 |
259 |
self.p_norm = p_norm |
self.p_norm = p_norm |
249 |
|
|
|
250 |
260 |
self.ns = ns |
self.ns = ns |
|
261 |
|
|
251 |
262 |
self.auto_update = auto_update |
self.auto_update = auto_update |
252 |
|
self.sample_callback = sample_callback |
|
|
263 |
|
self.on_add_mixed = on_add_mixed |
|
264 |
|
self.on_update_mixed = on_update_mixed |
|
265 |
|
self.on_delete_mixed = on_delete_mixed |
253 |
266 |
|
|
254 |
267 |
#č chcu, aby behem dalších iterací pulka (p_base) tečiček dopadala do buňky |
#č chcu, aby behem dalších iterací pulka (p_base) tečiček dopadala do buňky |
255 |
268 |
self.base_r = sball.Sball(sample_box.nvar).get_r(1 - p_base) |
self.base_r = sball.Sball(sample_box.nvar).get_r(1 - p_base) |
256 |
269 |
|
|
|
270 |
|
|
257 |
271 |
self._nsim = 0 |
self._nsim = 0 |
258 |
|
self.nodes = {} |
|
|
272 |
|
self.nodes = WardrobeDict() |
|
273 |
|
# checklist |
|
274 |
|
self.couples = {} #č |
259 |
275 |
self.direct_contacts = {} |
self.direct_contacts = {} |
260 |
276 |
self.couples_stats = {} |
self.couples_stats = {} |
261 |
277 |
self.mixed_couples = set() |
self.mixed_couples = set() |
262 |
278 |
self.red_couples = set() |
self.red_couples = set() |
263 |
279 |
#č uvnitř třídy ještě můžeme ordnung dodržovat |
#č uvnitř třídy ještě můžeme ordnung dodržovat |
264 |
280 |
self._indices_to_update = set() |
self._indices_to_update = set() |
265 |
|
self._couples_to_check_outsides = [] |
|
266 |
|
|
|
|
281 |
|
self._nodes_to_check_outsides = [] |
267 |
282 |
|
|
268 |
283 |
|
|
269 |
284 |
self.update() |
self.update() |
|
285 |
|
|
|
286 |
|
|
270 |
287 |
|
|
|
288 |
|
""" |
|
289 |
|
## checklist |
|
290 |
|
self.couples = {} # dict of active contacts {couple_indices: list of nodes} |
|
291 |
|
self.direct_contacts = {} # dict of direct contacts {couple_indices: node_idx} |
|
292 |
|
self.couples_stats = {} |
|
293 |
|
self.mixed_couples = set() |
|
294 |
|
self.red_couples = set() |
|
295 |
|
self._add_indices_to_update(j) # set of all the "j"s for the mutual updates |
|
296 |
|
self._nodes_to_check_outsides = [] #č seznam indexů, projedu jednoduše smyčkou |
|
297 |
|
""" |
271 |
298 |
|
|
272 |
299 |
def _update(self): |
def _update(self): |
273 |
300 |
if self.auto_update: |
if self.auto_update: |
|
... |
... |
class ContactVoronoi: |
291 |
318 |
|
|
292 |
319 |
self.sampled_plan_model = getattr(self.sample_box, self.model_space) |
self.sampled_plan_model = getattr(self.sample_box, self.model_space) |
293 |
320 |
self.tree = KDTree(sampled_plan_model) |
self.tree = KDTree(sampled_plan_model) |
294 |
|
|
|
|
321 |
|
self.CS = ConvexSolver(sampled_plan_model) |
295 |
322 |
|
|
296 |
323 |
|
|
297 |
324 |
#č já vím, že sample box je pokážde přepočítavá |
#č já vím, že sample box je pokážde přepočítavá |
|
... |
... |
class ContactVoronoi: |
303 |
330 |
#č a omezujeme se pouse nejbližšími vzorky |
#č a omezujeme se pouse nejbližšími vzorky |
304 |
331 |
#č disjunktnost je pořad zajištěna |
#č disjunktnost je pořad zajištěna |
305 |
332 |
#č a môžeme všechny nasbírané pravděpodobnosti jednoduše sčítat |
#č a môžeme všechny nasbírané pravděpodobnosti jednoduše sčítat |
|
333 |
|
# good old event ids... |
|
334 |
|
# -1 = 'outside', 0=success, 1=failure, 2=mix |
306 |
335 |
for i in range(self._nsim, nsim): |
for i in range(self._nsim, nsim): |
307 |
336 |
if i in failure_points: #č první je červený |
if i in failure_points: #č první je červený |
308 |
337 |
for j in range(i): |
for j in range(i): |
309 |
338 |
if j in failure_points: |
if j in failure_points: |
310 |
339 |
#č červený kontakt |
#č červený kontakt |
311 |
|
self.red_couples.add((i, j)) |
|
312 |
|
self.onboard_couple((i, j)) |
|
|
340 |
|
# -1 = 'outside', 0=success, 1=failure, 2=mix |
|
341 |
|
self.onboard_couple((i, j), event_id=1) |
313 |
342 |
else: |
else: |
314 |
343 |
#č žlutý kontakt |
#č žlutý kontakt |
315 |
|
self.mixed_couples.add((i, j)) |
|
316 |
|
self.onboard_couple((i, j)) |
|
|
344 |
|
# -1 = 'outside', 0=success, 1=failure, 2=mix |
|
345 |
|
self.onboard_couple((i, j), event_id=2) |
317 |
346 |
|
|
318 |
347 |
else: #č první je zelený |
else: #č první je zelený |
319 |
348 |
for j in range(i): |
for j in range(i): |
320 |
349 |
if j in failure_points: |
if j in failure_points: |
321 |
350 |
#č žlutý kontakt |
#č žlutý kontakt |
322 |
|
self.mixed_couples.add((i, j)) |
|
323 |
|
self.onboard_couple((i, j)) |
|
|
351 |
|
# -1 = 'outside', 0=success, 1=failure, 2=mix |
|
352 |
|
self.onboard_couple((i, j), event_id=2) |
324 |
353 |
|
|
325 |
354 |
#else: #č jinak nič |
#else: #č jinak nič |
326 |
355 |
#оӵ Вож. Кулэ ӧвӧл |
#оӵ Вож. Кулэ ӧвӧл |
|
... |
... |
class ContactVoronoi: |
346 |
375 |
|
|
347 |
376 |
self._nsim = nsim |
self._nsim = nsim |
348 |
377 |
|
|
|
378 |
|
def _add_indices_to_update(self, j): |
|
379 |
|
#č i-té indexy jsou ty čerstvé, přes ně iterujeme |
|
380 |
|
#č j-té - ty bežné, protí ním rozhodujeme |
|
381 |
|
if j < self._nsim: |
|
382 |
|
self._indices_to_update.add(j) |
|
383 |
|
|
349 |
384 |
|
|
350 |
|
def update_couple(self, couple_indices): |
|
|
385 |
|
def update_couple(self, couple_indices, event_id): |
|
386 |
|
pass |
|
387 |
|
|
|
388 |
|
def invalidate_couple(self, couple_indices): |
351 |
389 |
pass |
pass |
352 |
390 |
|
|
|
391 |
|
def onboard_couple(self, couple_indices, event_id): |
|
392 |
|
""" |
|
393 |
|
č Obecně, skoro pro každou dimenzi v rozmězí od 2D do 100D |
|
394 |
|
č nejlevnější bylo |
|
395 |
|
č 1. strom s pár bodíků |
|
396 |
|
č 2. ConvexSolver |
|
397 |
|
č 3. strom s nsim bodíků |
|
398 |
|
""" |
353 |
399 |
|
|
354 |
|
def onboard_couple(self, couple_indices): |
|
355 |
400 |
#č (i, j)? poprvé ta čísla vidíme |
#č (i, j)? poprvé ta čísla vidíme |
356 |
401 |
i, j = couple_indices |
i, j = couple_indices |
357 |
402 |
|
|
|
403 |
|
# 1. Nejdřív - zda je přímy kontakt, prostě na usečce |
358 |
404 |
contact_node_model = np.sum(self.sampled_plan_model[[i, j]], axis=0) |
contact_node_model = np.sum(self.sampled_plan_model[[i, j]], axis=0) |
359 |
405 |
np.divide(contact_node_model, 2, out=contact_node_model) |
np.divide(contact_node_model, 2, out=contact_node_model) |
360 |
406 |
|
|
|
407 |
|
dd, ii = tree.query(contact_node_model, k=2, p=self.p_norm) |
|
408 |
|
#č taková kontola je rychlejší jak volání numpy množinových funkcí |
|
409 |
|
if (i in ii) and (j in ii): |
|
410 |
|
|
|
411 |
|
#č přídat, zavolat kólbek |
|
412 |
|
f = self.sample_box.f_model |
|
413 |
|
#č zde je jistá nadbytečnost. |
|
414 |
|
#č Pro třídu je jednodušší prostě pokažde spočítat bodík, |
|
415 |
|
#č než je ukladat, sledovat, invalidovat. |
|
416 |
|
#č Jedině to může být zajimavé pro skřiňku |
|
417 |
|
#č a i pro ně - jenom smíšené bodíky |
|
418 |
|
contact_node = f.new_sample(contact_node_model, self.model_space) |
|
419 |
|
contact_node = CandyBox(contact_node, \ |
|
420 |
|
dd1=dd[:, 0], dd2=dd[:, 1], ii1=ii[:, 0], ii2=ii[:, 1]) |
|
421 |
|
|
|
422 |
|
idx = self.nodes.add(contact_node) |
|
423 |
|
## checklist |
|
424 |
|
#self.couples #č další kód, kdyby něco |
|
425 |
|
self.direct_contacts[couple_indices] = idx |
|
426 |
|
#self._nodes_to_check_outsides #č nemůže být outside |
|
427 |
|
|
|
428 |
|
# -1 = 'outside', 0=success, 1=failure, 2=mix |
|
429 |
|
if (event_id == 2) and (self.on_add_mixed is not None): |
|
430 |
|
self.on_add_mixed(idx) |
|
431 |
|
|
|
432 |
|
|
|
433 |
|
|
|
434 |
|
#č přímy kontakt není, ale snad i tak existuje? |
|
435 |
|
elif not self.CS.is_couple(couple_indices): |
|
436 |
|
#č neexistuje, netřeba nic robiť, j-ko i-čkem není dotčen |
|
437 |
|
return |
|
438 |
|
|
|
439 |
|
#č kontakt existuje, deme vzorkovat |
|
440 |
|
|
|
441 |
|
|
|
442 |
|
|
|
443 |
|
# checklist |
|
444 |
|
#self.couples = #č přídá sample_couple, když se mu povede si oblast chytit |
|
445 |
|
#self.direct_contacts #č není |
|
446 |
|
#self.couples_stats = {} |
|
447 |
|
#self.mixed_couples = set() |
|
448 |
|
#self.red_couples = set() |
|
449 |
|
self._add_indices_to_update(j) #č kontakt je, apdejtneme |
|
450 |
|
#self._nodes_to_check_outsides #č přídá sample_couple kdyby něco |
|
451 |
|
|
|
452 |
|
def IS_norm(f, mean=0, std=1, sampling_space='G', nis=1000, design=None): |
|
453 |
|
""" |
|
454 |
|
mean: [0.05, 2, 100500] |
|
455 |
|
std: [0.05, 2, 100500] |
|
456 |
|
|
|
457 |
|
design(nis, nvar) should return sampling plan in hypercube (U space)! |
|
458 |
|
""" |
|
459 |
|
|
|
460 |
|
sampling_plan_N, pdf = IS_stat.get_norm_plan(nis, f.nvar, mean, std, design) |
|
461 |
|
|
|
462 |
|
#č tady musíme provést jeden trik |
|
463 |
|
#č totež jako v IS_like - ve výsledku dycky dostaneme f_model |
|
464 |
|
to_sample = f.new_sample(sampling_plan_N, sampling_space) #č naše N-ko smerdžíme ako G-čko |
|
465 |
|
w = to_sample.pdf(sampling_space) / pdf #č snad je to správně |
|
466 |
|
|
|
467 |
|
#č vahy máme |
|
468 |
|
#č zabalme do boxu |
|
469 |
|
#č zbytek už nejsou naši starosti |
|
470 |
|
return CandyBox(to_sample, w=w) |
|
471 |
|
|
|
472 |
|
|
|
473 |
|
|
|
474 |
|
def sample_couple(self, couple_indices, event_id, center, r): |
|
475 |
|
"""č máme nový manželský páreček |
|
476 |
|
оӵ Выль кузо. Мар меда кароно? |
|
477 |
|
č zde jenom vzorkujeme; |
|
478 |
|
č integrování, vyhodnocování je jinde a později! |
|
479 |
|
""" |
|
480 |
|
|
|
481 |
|
|
|
482 |
|
|
|
483 |
|
# use IS sampling density with center equal to the current point |
|
484 |
|
# set the minimum distance as the standard deviation of IS densisty |
|
485 |
|
# mindist dělíme dvěma - hranice buňky je na půlcesty |
|
486 |
|
# a dělíme base_r'kem k získání sigmy. Přemejmenším empiricky mně to vychází |
|
487 |
|
# (u Norm'u vždy zadaváme směrodatnou odchylku) |
|
488 |
|
h_i = [stats.norm(sampled_plan_sing[i,j], mindist_sing[i]/2/base_r) for j in range(sample_box.nvar)] |
|
489 |
|
h = f_models.UnCorD(h_i) |
|
490 |
|
|
|
491 |
|
|
|
492 |
|
# select nis points from IS density |
|
493 |
|
# sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku |
|
494 |
|
h_plan = IS(f, h, space_from_h='R', space_to_f=sampling_space, Nsim=nis) |
|
495 |
|
|
|
496 |
|
# součet váh nemá cenu kontrolovat, jedná tam nebude, nevyjde |
|
497 |
|
|
|
498 |
|
h_plan_model = getattr(h_plan, model_space) |
|
499 |
|
dd, ii = tree.query(h_plan_model, k=1, p=p_norm) |
|
500 |
|
|
|
501 |
|
|
|
502 |
|
# nechám s velkým písmenkem |
|
503 |
|
Vor_mask = ii==i |
|
504 |
|
h_plan_model_ma = h_plan_model[Vor_mask] |
|
505 |
|
weights_sim = h_plan.w[Vor_mask] |
|
506 |
|
# dd1 jsou vzdalenosti tečiček do centra Voroneho buňky |
|
507 |
|
dd1 = dd[Vor_mask] |
|
508 |
|
|
|
509 |
|
trace = 0 |
|
510 |
|
nis_ma = len(weights_sim) |
|
511 |
|
|
|
512 |
|
h_plan_sing_ma = getattr(h_plan, sampling_space)[Vor_mask] |
|
513 |
|
S_bc = np.cov(h_plan_sing_ma, rowvar=False, bias=True, aweights=weights_sim) |
|
514 |
|
barycenter = np.average(h_plan_sing_ma, axis=0, weights=weights_sim) |
|
515 |
|
# effective nis |
|
516 |
|
nis_eff = nis |
|
517 |
|
#print("start", i) |
|
518 |
|
|
|
519 |
|
while trace < np.trace(S_bc): |
|
520 |
|
# number of nodes inside the cell matters too |
|
521 |
|
trace = np.trace(S_bc)# * nis_ma |
|
522 |
|
#print(S_bc) |
|
523 |
|
#print(nis_ma) |
|
524 |
|
|
|
525 |
|
|
|
526 |
|
|
|
527 |
|
|
|
528 |
|
# matika |
|
529 |
|
w, v = np.linalg.eig(S_bc) |
|
530 |
|
|
|
531 |
|
# use IS sampling density with center equal to the cell's barycenter |
|
532 |
|
# set the minimum distance as the standard deviation of the IS densisty |
|
533 |
|
# u stats.norm zadáváme směrodatnou odchylku, to je asi správné |
|
534 |
|
sigmas = np.sqrt(w) / base_r #(sample_box.nvar+2) #! dosadit standard deviation podle chutí |
|
535 |
|
h_i = [stats.norm(0, sigmas[j]) for j in range(sample_box.nvar)] |
|
536 |
|
# rozdělení ve vlastním prostoru |
|
537 |
|
# select nis = 100 points from IS density |
|
538 |
|
h_L = f_models.UnCorD(h_i)(nis) |
|
539 |
|
|
|
540 |
|
# здесь уже так легко не отделаемся. Трансформовать кароно. |
|
541 |
|
h_plan_bc = (v @ h_L.R.T).T |
|
542 |
|
h_plan_sing = h_plan_bc + barycenter |
|
543 |
|
|
|
544 |
|
# sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku |
|
545 |
|
h_plan_part = f.new_sample(h_plan_sing, sampling_space) |
|
546 |
|
|
|
547 |
|
|
|
548 |
|
# jdeme na ty hustoty |
|
549 |
|
# mně příjde, že je to legalní |
|
550 |
|
# sice samply podporujou maskovaní, to je ale drahé |
|
551 |
|
weights_sim_part = h_plan_part.pdf(sampling_space) / h_L.pdf('R') # snad je to správně |
|
552 |
|
h_plan.add_sample(CandyBox(h_plan_part, w=weights_sim_part)) |
|
553 |
|
|
|
554 |
|
# vyfiltrujeme vzorky |
|
555 |
|
h_plan_model_part = getattr(h_plan_part, model_space) |
|
556 |
|
dd, ii = tree.query(h_plan_model_part, k=1, p=p_norm) |
|
557 |
|
|
|
558 |
|
# parta |
|
559 |
|
Vor_mask_part = ii==i |
|
560 |
|
weights_sim_part = weights_sim_part[Vor_mask_part] |
|
561 |
|
|
|
562 |
|
nis_ma = len(weights_sim_part) |
|
563 |
|
|
|
564 |
|
# zajišťovat Vor_mask je docela zbytečně, je to jen pro out_nodes, |
|
565 |
|
# které se zatím nikdě nepouživá |
|
566 |
|
Vor_mask = np.append(Vor_mask, Vor_mask_part) |
|
567 |
|
|
|
568 |
|
|
|
569 |
|
|
|
570 |
|
|
|
571 |
|
h_plan_model_ma = np.vstack((h_plan_model_ma, h_plan_model_part[Vor_mask_part])) |
|
572 |
|
weights_sim = np.append(weights_sim, weights_sim_part) |
|
573 |
|
# dd1 jsou vzdalenosti tečiček do centra Voroneho buňky |
|
574 |
|
dd1 = np.append(dd1, dd[Vor_mask_part]) |
|
575 |
|
|
|
576 |
|
# zkusím těžiště počitat jen pro partu - možná tak algoritmus bude agresivnější? |
|
577 |
|
#barycenter = np.average(h_plan_sing[Vor_mask_part], axis=0, weights=weights_sim_part) |
|
578 |
|
h_plan_sing_ma = np.vstack((h_plan_sing_ma, h_plan_sing[Vor_mask_part])) |
|
579 |
|
#S_bc = np.cov(h_plan_sing[Vor_mask_part], rowvar=False, bias=True, aweights=weights_sim_part) |
|
580 |
|
S_bc = np.cov(h_plan_sing_ma, rowvar=False, bias=True, aweights=weights_sim) |
|
581 |
|
barycenter = np.average(h_plan_sing_ma, axis=0, weights=weights_sim) |
|
582 |
|
nis_eff += nis |
|
583 |
|
|
|
584 |
|
|
|
585 |
|
|
|
586 |
|
#print(S_bc) |
|
587 |
|
#print(nis_ma) |
|
588 |
|
|
|
589 |
|
|
|
590 |
|
|
|
591 |
|
|
|
592 |
|
cell_stats = dict() |
|
593 |
|
# musí sa (na konci) rovnat jedne |
|
594 |
|
# opravdu dělíme nis'em, jako v normálním IS |
|
595 |
|
# nikoliv počtem příjatých bodíků, |
|
596 |
|
# protože průměrná vaha je o hodně mene významná metrika |
|
597 |
|
cell_stats['cell_probability'] = np.sum(weights_sim) / nis_eff |
|
598 |
|
|
361 |
599 |
|
|
362 |
|
# self.direct_contacts[(i, j)] = None |
|
363 |
|
# self.nodes[(i, j)] = None |
|
|
600 |
|
# tu bacha! |
|
601 |
|
# takhle se počíta, pokud se netrapíme gradijentem |
|
602 |
|
# a je to trošiňku optimizovaný, takže čert se nevyzná |
|
603 |
|
if gradient is None: |
|
604 |
|
# indexy ii nás moc nezajimajou |
|
605 |
|
# vzdalenosti snad byjsme zvladli použit? |
|
606 |
|
|
|
607 |
|
|
|
608 |
|
dd2, ii2 = tree.query(h_plan_model_ma, k=[2], p=p_norm) |
|
609 |
|
dd2 = dd2.flatten() |
|
610 |
|
ii2 = ii2.flatten() |
|
611 |
|
|
|
612 |
|
# tahle hračka s indexy je pro numpy poměrně drahá |
|
613 |
|
failsii_2 = failsi[ii2] |
|
614 |
|
# jeden vzorek (včetně hustoty PDF[i]) je nám vždy znám |
|
615 |
|
# porucha |
|
616 |
|
if failsi[i]: |
|
617 |
|
points_1 = PDF[i] * dd2 |
|
618 |
|
node_pf_estimations = points_1 / (points_1 + PDF[ii2] * dd1) |
|
619 |
|
node_pf_estimations = np.where(failsii_2,1, node_pf_estimations) |
|
620 |
|
node_pf_pure_estimations = dd2 / (dd1 + dd2) |
|
621 |
|
node_pf_pure_estimations = np.where(failsii_2,1, node_pf_pure_estimations) |
|
622 |
|
|
|
623 |
|
cell_stats['Voronoi_2_point_upper_bound'] = cell_stats['cell_probability'] |
|
624 |
|
cell_stats['Voronoi_2_point_failure_rate'] = np.sum(weights_sim * node_pf_estimations) / nis_eff |
|
625 |
|
cell_stats['Voronoi_2_point_pure_failure_rate'] = np.sum(weights_sim * node_pf_pure_estimations) / nis_eff |
|
626 |
|
cell_stats['Voronoi_2_point_lower_bound'] = np.sum(weights_sim[failsii_2]) / nis_eff |
|
627 |
|
cell_stats['Voronoi_failure_rate'] = cell_stats['cell_probability'] |
|
628 |
|
nodes=CandyBox(h_plan.sampling_plan[Vor_mask], w=weights_sim, node_pf_estimations=node_pf_estimations,\ |
|
629 |
|
node_pf_pure_estimations=node_pf_pure_estimations, dd1=dd1, dd2=dd2, ii2=ii2) |
|
630 |
|
|
|
631 |
|
# neporucha |
|
632 |
|
else: |
|
633 |
|
dd1 = dd1[failsii_2] |
|
634 |
|
dd2 = dd2[failsii_2] |
|
635 |
|
points_1 = PDF[i] * dd2 |
|
636 |
|
points_2 = PDF[ii2[failsii_2]] * dd1 |
|
637 |
|
|
|
638 |
|
node_pf_estimations = points_2 / (points_1 + points_2) |
|
639 |
|
node_pf_pure_estimations = dd1 / (dd1 + dd2) |
|
640 |
|
|
|
641 |
|
cell_stats['Voronoi_2_point_upper_bound'] = np.sum(weights_sim[failsii_2]) / nis_eff |
|
642 |
|
cell_stats['Voronoi_2_point_failure_rate'] = np.sum(weights_sim[failsii_2]*node_pf_estimations) / nis_eff |
|
643 |
|
cell_stats['Voronoi_2_point_pure_failure_rate'] = np.sum(weights_sim[failsii_2] * node_pf_pure_estimations) / nis_eff |
|
644 |
|
cell_stats['Voronoi_2_point_lower_bound'] = 0 |
|
645 |
|
cell_stats['Voronoi_failure_rate'] = 0 |
|
646 |
|
nodes=CandyBox(h_plan.sampling_plan[Vor_mask][failsii_2], w=weights_sim[failsii_2], node_pf_estimations=node_pf_estimations,\ |
|
647 |
|
node_pf_pure_estimations=node_pf_pure_estimations, dd1=dd1, dd2=dd2, ii2=ii2[failsii_2]) |
|
648 |
|
|
|
649 |
|
# take something with corresponding length |
|
650 |
|
zeros = np.zeros(len(weights_sim) - len(dd2)) |
|
651 |
|
# add remaining nodes |
|
652 |
|
nodes.add_sample(CandyBox(h_plan.sampling_plan[Vor_mask][~failsii_2], w=weights_sim[~failsii_2], node_pf_estimations=zeros,\ |
|
653 |
|
node_pf_pure_estimations=zeros, ii2=ii2[~failsii_2])) |
364 |
654 |
|
|
365 |
655 |
|
|
366 |
|
dd, ii = tree.query(h_plan_model_part, k=1, p=p_norm) |
|
|
656 |
|
# takhle - pokud chceme gradient použit |
|
657 |
|
# je třeba eště zoptimalizovať |
|
658 |
|
else: |
|
659 |
|
# kolik bodíků jsou nejbližší k mému vzorečkovi |
|
660 |
|
# np.empty() implicitně má dtype=float |
|
661 |
|
# tyhle blbosti ponechám jen kvůli callbackovi |
|
662 |
|
node_pf_estimations = np.empty(len(h_plan_model_ma)) |
|
663 |
|
node_pf_pure_estimations = np.empty(len(h_plan_model_ma))# pure distance estimation |
|
664 |
|
node_failsi = np.empty(len(h_plan_model_ma), dtype=np.bool) # for L1 Voronoi # co to je za L1 Voronoi? |
|
665 |
|
|
|
666 |
|
# projdeme přes každej bodíček |
|
667 |
|
for node_idx in range(len(h_plan_model_ma)): |
|
668 |
|
# KDTree byl použit jen k rozdělení na disjunktní úseky, veškerej děj se odehravá tu |
|
669 |
|
# a to je všechno kvůli gradientu |
|
670 |
|
node = h_plan_model_ma[node_idx] |
|
671 |
|
# axis=1 - sčítá všechy směry dohromady, vysledkem je 1D pole rozměru nsim |
|
672 |
|
inode2points_model_matrix = np.sum(np.abs(((sampled_plan_model - node) * gradient(node))**p_norm), axis=1) |
|
673 |
|
#print(inode2points_Rd_matrix) |
|
674 |
|
|
|
675 |
|
""" |
|
676 |
|
partition - |
|
677 |
|
Creates a copy of the array with its elements rearranged in such a way that |
|
678 |
|
the value of the element in k-th position is in the position it would be in a sorted array. |
|
679 |
|
All elements smaller than the k-th element are moved before this element |
|
680 |
|
and all equal or greater are moved behind it. The ordering of the elements in the two partitions is undefined. |
|
681 |
|
""" |
|
682 |
|
idx = np.argpartition(inode2points_model_matrix, 1) # musí tu bejt 1, počítá sa od nuly |
|
683 |
|
# je to správný, neboť numpy zaručuje, že druhej prvek (s indexem 1) bude na druhem místě |
|
684 |
|
node_failsi[node_idx] = failsi[idx[0]] |
|
685 |
|
|
|
686 |
|
|
|
687 |
|
points_weight = PDF[idx[0]] / inode2points_model_matrix[idx[0]] + PDF[idx[1]] / inode2points_model_matrix[idx[1]] |
|
688 |
|
points_distances = 1 / inode2points_model_matrix[idx[0]] + 1 / inode2points_model_matrix[idx[1]] |
|
689 |
|
|
|
690 |
|
failure_weight = int(failsi[idx[0]]) * PDF[idx[0]] / inode2points_model_matrix[idx[0]] |
|
691 |
|
failure_weight += int(failsi[idx[1]]) * PDF[idx[1]] / inode2points_model_matrix[idx[1]] |
367 |
692 |
|
|
|
693 |
|
failure_distance = int(failsi[idx[0]]) / inode2points_model_matrix[idx[0]] + int(failsi[idx[1]]) / inode2points_model_matrix[idx[1]] |
|
694 |
|
|
|
695 |
|
|
|
696 |
|
node_pf_estimations[node_idx] = failure_weight/points_weight |
|
697 |
|
node_pf_pure_estimations[node_idx] = failure_distance/points_distances |
|
698 |
|
|
|
699 |
|
|
|
700 |
|
cell_stats['Voronoi_2_point_upper_bound'] = np.sum(h_plan.w[Vor_mask]*np.ceil(node_pf_estimations)) / nis_eff |
|
701 |
|
cell_stats['Voronoi_2_point_failure_rate'] = np.sum(h_plan.w[Vor_mask]*node_pf_estimations) / nis_eff |
|
702 |
|
cell_stats['Voronoi_2_point_pure_failure_rate'] = np.sum(h_plan.w[Vor_mask]*node_pf_pure_estimations) / nis_eff |
|
703 |
|
cell_stats['Voronoi_2_point_lower_bound'] = np.sum(h_plan.w[Vor_mask]*np.floor(node_pf_estimations)) / nis_eff |
|
704 |
|
cell_stats['Voronoi_failure_rate'] = np.sum(h_plan.w[Vor_mask]*node_failsi) / nis_eff |
|
705 |
|
|
|
706 |
|
nodes=CandyBox(h_plan.sampling_plan[Vor_mask], w=h_plan.w[Vor_mask], node_pf_estimations=node_pf_estimations,\ |
|
707 |
|
node_pf_pure_estimations=node_pf_pure_estimations, node_failsi=node_failsi) |
|
708 |
|
|
|
709 |
|
|
|
710 |
|
|
|
711 |
|
for key, value in cell_stats.items(): |
|
712 |
|
global_stats[key] += value |
|
713 |
|
|
|
714 |
|
# kolbek ↓ |
|
715 |
|
if failsi[i]: |
|
716 |
|
cell_stats['event'] = 'failure' |
|
717 |
|
else: |
|
718 |
|
cell_stats['event'] = 'success' |
|
719 |
|
callback(estimation=estimation, nodes=nodes, cell_stats=cell_stats, out_nodes=h_plan[~Vor_mask]) |
|
720 |
|
|
|
721 |
|
|
|
722 |
|
|
|
723 |
|
|
368 |
724 |
# find distance to the nearest sampling point (from all points) |
# find distance to the nearest sampling point (from all points) |
369 |
725 |
dd2, ii2 = tree_sampling.query(sampled_plan_sing, k=[2], p=2) |
dd2, ii2 = tree_sampling.query(sampled_plan_sing, k=[2], p=2) |
370 |
726 |
mindist_sing = dd2.flatten() |
mindist_sing = dd2.flatten() |
371 |
727 |
|
|
372 |
728 |
|
|
|
729 |
|
def _store_masked(self, nodes_model, couple_indices, event_id): |
|
730 |
|
""" |
|
731 |
|
č logika celé třídy je, že máme-li příslušné kontaktu bodíky, |
|
732 |
|
č tak stojí za to je uložit. |
|
733 |
|
č Proto metod maskuje a je-li co ukladat tak i rovnou uloži |
|
734 |
|
č a vratí jejich index (jinak vratí nulu). |
|
735 |
|
č To ale vyžaduje trochu větší pozornost - |
|
736 |
|
č bodíky pak někdo má odebrat a invalidovat. |
|
737 |
|
č Proto i podtržitko v názvu |
|
738 |
|
""" |
|
739 |
|
i, j = couple_indices |
|
740 |
|
dd, ii = tree.query(nodes_model, k=2, p=self.p_norm) |
|
741 |
|
#č pro jediný bod jednoducha kontrola je rychlejší |
|
742 |
|
#č jak volání numpy množinových funkcí |
|
743 |
|
#if (i in ii) and (j in ii): |
|
744 |
|
|
|
745 |
|
# similar approach is used in numpy code |
|
746 |
|
mask = np.all((ii == i) | (ii == j), axis=1) |
|
747 |
|
|
|
748 |
|
if not np.any(mask): |
|
749 |
|
return False |
373 |
750 |
|
|
374 |
751 |
|
|
375 |
|
#č není znamená, že nikdy ani nebylo |
|
376 |
|
# if ((i, j) in self.nodes) or ((i, j) in self.direct_contacts): |
|
377 |
|
# #č None znamená, že kontakt je pro nás pase |
|
378 |
|
# if self.nodes[(i, j)] is not None: |
|
379 |
|
# self.update_couple((i, j)) |
|
380 |
|
# |
|
381 |
|
# else: #č (i, j)? poprvé ta čísla vidíme |
|
382 |
|
# self.onboard_couple((i, j)) |
|
|
752 |
|
#č přídat, zavolat kólbek |
|
753 |
|
f = self.sample_box.f_model |
|
754 |
|
nodes = f.new_sample(nodes_model[mask], self.model_space) |
|
755 |
|
nodes = CandyBox(nodes, dd1=dd[:, 0], dd2=dd[:, 1], \ |
|
756 |
|
ii1=ii[:, 0], ii2=ii[:, 1]) |
383 |
757 |
|
|
|
758 |
|
idx = self.nodes.add(nodes) |
|
759 |
|
## checklist |
|
760 |
|
#self.couples #č další kód, kdyby něco |
|
761 |
|
self.direct_contacts[couple_indices] = idx |
|
762 |
|
#self._nodes_to_check_outsides #č nemůže být outside |
|
763 |
|
|
|
764 |
|
# -1 = 'outside', 0=success, 1=failure, 2=mix |
|
765 |
|
if (event_id == 2) and (self.on_add_mixed is not None): |
|
766 |
|
self.on_add_mixed(idx) |
|
767 |
|
|
|
768 |
|
## checklist |
|
769 |
|
self.couples = {} # dict of active contacts {couple_indices: list of nodes} |
|
770 |
|
self.direct_contacts = {} # dict of direct contacts {couple_indices: node_idx} |
|
771 |
|
self.couples_stats = {} |
|
772 |
|
self.mixed_couples = set() |
|
773 |
|
self.red_couples = set() |
|
774 |
|
self._add_indices_to_update(j) # set of all the "j"s for the mutual updates |
|
775 |
|
self._nodes_to_check_outsides = [] #č seznam indexů, projedu jednoduše smyčkou |
|
776 |
|
|
|
777 |
|
self.red_couples.add((i, j)) |
|
778 |
|
|
|
779 |
|
self.mixed_couples.add((i, j)) |
384 |
780 |
|
|
385 |
|
#č i-té indexy jsou ty čerstvé, přes ně iterujeme |
|
386 |
|
#č j-té - ty bežné, protí ním rozhodujeme |
|
387 |
|
if j < self._nsim: |
|
388 |
|
self._indices_to_update.add(j) |
|
389 |
781 |
|
|
390 |
782 |
|
|
391 |
783 |
|
|
|
784 |
|
|
|
785 |
|
|
|
786 |
|
|
|
787 |
|
|
392 |
788 |
def integrate_couple(self, couple_indices): |
def integrate_couple(self, couple_indices): |
393 |
789 |
"""č máme nový manželský páreček |
"""č máme nový manželský páreček |
394 |
790 |
оӵ Выль кузо. Мар меда кароно? |
оӵ Выль кузо. Мар меда кароно? |
|
... |
... |
class ContactVoronoi: |
632 |
1028 |
cell_stats['event'] = 'failure' |
cell_stats['event'] = 'failure' |
633 |
1029 |
else: |
else: |
634 |
1030 |
cell_stats['event'] = 'success' |
cell_stats['event'] = 'success' |
635 |
|
callback(estimation=estimation, nodes=nodes, cell_stats=cell_stats, out_nodes=h_plan[~Vor_mask]) |
|
636 |
|
|
|
637 |
|
|
|
638 |
|
|
|
639 |
|
|
|
|
1031 |
|
callback(estimation=estimation, nodes=nodes, cell_stats=cell_stats, out_nodes=h_plan[~Vor_mask]) |
640 |
1032 |
|
|
641 |
1033 |
|
|
642 |
1034 |
|
|