File wellmet/voronoi.py changed (mode: 100644) (index 00974a5..b0f3508) |
... |
... |
from .IS_stat import IS |
20 |
20 |
from .candynodes import CandyNodes |
from .candynodes import CandyNodes |
21 |
21 |
from . import sball |
from . import sball |
22 |
22 |
|
|
|
23 |
|
##č (jednotka dočasnosti - jeden furt) |
23 |
24 |
|
|
24 |
25 |
|
|
25 |
26 |
|
|
|
... |
... |
class WardrobeDict(dict): |
40 |
41 |
|
|
41 |
42 |
|
|
42 |
43 |
TData = collections.namedtuple('TData', ('X', 'imask', 'dd')) |
TData = collections.namedtuple('TData', ('X', 'imask', 'dd')) |
|
44 |
|
#VorSTM = collections.namedtuple('VorSTM', \ |
|
45 |
|
# ('nsim', 'failure', 'mixed', 'pfv', 'pfw')) |
|
46 |
|
|
|
47 |
|
## TRI-compatible estimation |
|
48 |
|
#TriSTM = collections.namedtuple('TriSTM', \ |
|
49 |
|
# ('outside', 'success', 'failure', 'mixed', 'pfv', 'pfw')) |
|
50 |
|
|
|
51 |
|
|
43 |
52 |
|
|
44 |
53 |
# HalfspaceVoronoi |
# HalfspaceVoronoi |
45 |
54 |
# TwoFaces :) |
# TwoFaces :) |
|
... |
... |
class ContactVoronoi: |
595 |
604 |
#č já chcu, aby třída byla blbuvzdorná |
#č já chcu, aby třída byla blbuvzdorná |
596 |
605 |
#č aby šlo ji vytvořit na prazdné skřiňce |
#č aby šlo ji vytvořit na prazdné skřiňce |
597 |
606 |
#č aby neotravovala chybáma |
#č aby neotravovala chybáma |
598 |
|
self._update() |
|
|
607 |
|
#self._update() |
599 |
608 |
|
|
600 |
609 |
|
|
601 |
610 |
# dump methods |
# dump methods |
|
... |
... |
class ContactVoronoi: |
616 |
625 |
self.couples = {} # dict of active contacts {couple_indices: list of nodes} |
self.couples = {} # dict of active contacts {couple_indices: list of nodes} |
617 |
626 |
self.couples_stats = {} |
self.couples_stats = {} |
618 |
627 |
self.mixed_couples = {} ##č pár slovníků reprezentativních sad vzorků |
self.mixed_couples = {} ##č pár slovníků reprezentativních sad vzorků |
619 |
|
self.red_couples = {} ## {couple_indices: node_idx} |
|
|
628 |
|
self.red_couples = {} ## {couple_indices: nodes} |
620 |
629 |
self._add_indices_to_update(j) |
self._add_indices_to_update(j) |
621 |
630 |
self._valid_outsides = set() #č zabyvá se tím bezprostředně integrovačka |
self._valid_outsides = set() #č zabyvá se tím bezprostředně integrovačka |
622 |
631 |
|
|
|
... |
... |
class ContactVoronoi: |
624 |
633 |
Nodes pipeline: |
Nodes pipeline: |
625 |
634 |
0. Sampling: generate coordinates |
0. Sampling: generate coordinates |
626 |
635 |
1. _store_masked(): Filter out, create CandyNodes |
1. _store_masked(): Filter out, create CandyNodes |
627 |
|
and (optionally) assign dd1, dd2, ii1, ii2 |
|
|
636 |
|
with d2 and (optionally) imask assigned |
628 |
637 |
2. assign weights (w) such as sum(w)=1 over entire domain |
2. assign weights (w) such as sum(w)=1 over entire domain |
629 |
|
3. score(): assign eig, barycenter and trace-based so called "score" |
|
|
638 |
|
3. score(): assign so called "score" |
630 |
639 |
""" |
""" |
631 |
640 |
|
|
|
641 |
|
|
|
642 |
|
#č funkce se nemá dělat toho víc, |
|
643 |
|
#č pokud by ji někdo opakovaně volal pořad dokola |
632 |
644 |
def _update(self): |
def _update(self): |
633 |
645 |
if self.auto_update: |
if self.auto_update: |
634 |
646 |
if self.basic_checks(): |
if self.basic_checks(): |
635 |
|
self.update_tree() |
|
636 |
647 |
##č můžeme udělat obraceně - |
##č můžeme udělat obraceně - |
637 |
648 |
##č je-li to nutné, nechť si uživatel třídy |
##č je-li to nutné, nechť si uživatel třídy |
638 |
649 |
##č zajístí správnost autsajdů |
##č zajístí správnost autsajdů |
|
... |
... |
class ContactVoronoi: |
648 |
659 |
return True |
return True |
649 |
660 |
elif self.sample_box.nsim < 2: |
elif self.sample_box.nsim < 2: |
650 |
661 |
return False |
return False |
651 |
|
elif not np.any(self.failsi): |
|
|
662 |
|
elif not np.any(self.sample_box.failsi): |
652 |
663 |
return False |
return False |
653 |
664 |
else: |
else: |
654 |
665 |
self._enabled = True |
self._enabled = True |
|
... |
... |
class ContactVoronoi: |
657 |
668 |
def invalidate_outside(self): |
def invalidate_outside(self): |
658 |
669 |
self._valid_outsides.clear() |
self._valid_outsides.clear() |
659 |
670 |
|
|
660 |
|
|
|
|
671 |
|
#č funkce se nemá dělat toho víc, |
|
672 |
|
#č pokud by ji někdo opakovaně volal pořad dokola |
661 |
673 |
def update_contacts(self): |
def update_contacts(self): |
662 |
674 |
#č já vím, že sample box pokážde failsi přepočítavá |
#č já vím, že sample box pokážde failsi přepočítavá |
663 |
675 |
self.failsi = failsi = self.sample_box.failsi |
self.failsi = failsi = self.sample_box.failsi |
664 |
|
|
|
665 |
676 |
self.points = getattr(self.sample_box, self.model_space) |
self.points = getattr(self.sample_box, self.model_space) |
666 |
677 |
|
|
667 |
|
self.CS = ConvexSpline(sampled_plan_model) |
|
668 |
|
self.PDF = self.sample_box.pdf(self.model_space) |
|
669 |
|
|
|
670 |
678 |
assert len(self._indices_to_update) == 0 |
assert len(self._indices_to_update) == 0 |
671 |
679 |
|
|
672 |
680 |
nis = self.ns |
nis = self.ns |
|
... |
... |
class ContactVoronoi: |
710 |
718 |
#č ale bez bodíků se aktualizace brzy skončí :) |
#č ale bez bodíků se aktualizace brzy skončí :) |
711 |
719 |
self.update_couple(couple_indices) |
self.update_couple(couple_indices) |
712 |
720 |
self._indices_to_update.clear() |
self._indices_to_update.clear() |
713 |
|
|
|
714 |
|
#č neplatí. Bohužel |
|
715 |
|
# #č pokud _nsim > 0, tak není možné, že by nebyly apdejty |
|
716 |
|
# #č "nejmenší" pár je (1, 0), takže |
|
717 |
|
# #č "i" může být v rozmězí od 1 do _nsim-1 (včetně) |
|
718 |
|
# #č u prvního indexu musíme dodržovat pořadí procházky |
|
719 |
|
# for i in range(self._nsim - 1, 0, -1): |
|
720 |
|
# if i in self._indices_to_update: |
|
721 |
|
# self._indices_to_update.remove(i) |
|
722 |
|
# for j in self._indices_to_update: |
|
723 |
|
# if ((i, j) in self.mixed_couples) or ((i, j) in self.red_couples): |
|
724 |
|
# self.update_couple((i, j)) |
|
725 |
|
# |
|
726 |
|
# #č nechcu ani kontrolovat zda-li, jestli, není-li. Mouchy. |
|
727 |
|
# #if |
|
728 |
|
# # break |
|
729 |
|
|
|
730 |
|
|
|
731 |
|
|
|
732 |
721 |
self._nsim = nsim |
self._nsim = nsim |
733 |
722 |
|
|
|
723 |
|
|
734 |
724 |
def _add_indices_to_update(self, j): |
def _add_indices_to_update(self, j): |
735 |
725 |
#č i-té indexy jsou ty čerstvé, přes ně iterujeme |
#č i-té indexy jsou ty čerstvé, přes ně iterujeme |
736 |
726 |
#č j-té - ty bežné, protí ním rozhodujeme |
#č j-té - ty bežné, protí ním rozhodujeme |
|
... |
... |
class ContactVoronoi: |
877 |
867 |
#č To ale nic neříká o hustotách |
#č To ale nic neříká o hustotách |
878 |
868 |
nodes.w = w * np.isfinite(w) |
nodes.w = w * np.isfinite(w) |
879 |
869 |
self.score(nodes) |
self.score(nodes) |
880 |
|
|
|
881 |
|
|
|
882 |
|
#č doporučíme (dočasně) k integraci |
|
883 |
|
self._recommend_to_integration(idx) |
|
|
870 |
|
|
|
871 |
|
|
|
872 |
|
|
884 |
873 |
|
|
885 |
874 |
#č jdeme do adaptivního IS vzorkování |
#č jdeme do adaptivního IS vzorkování |
886 |
875 |
#č zde jenom vzorkujeme; |
#č zde jenom vzorkujeme; |
887 |
876 |
#č integrování, vyhodnocování je jinde a později! |
#č integrování, vyhodnocování je jinde a později! |
888 |
|
|
|
889 |
|
nodes = self.sample_alike(tp, nodes, nis) |
|
|
877 |
|
score = -np.inf |
|
878 |
|
while nodes.score > score: |
|
879 |
|
score = nodes.score |
|
880 |
|
#č doporučíme (dočasně) k integraci |
|
881 |
|
self._recommend_to_integration(nodes) |
|
882 |
|
nodes = self.sample_alike(tp, nodes, nis) |
890 |
883 |
|
|
891 |
|
idx, dd, ii, mask = self._store_masked(nodes_model, couple_indices, event_id, h_pdf) |
|
892 |
|
if not idx: |
|
893 |
|
return None |
|
894 |
|
|
|
895 |
|
#č uspěch. Dodáme bodíkům vahy |
|
896 |
|
nodes = self.nodes[idx] |
|
897 |
|
w = nodes.pdf(self.model_space) / h_pdf[mask] / self.ns |
|
898 |
|
#č Tamta filtrace stromem zaručuje, |
|
899 |
|
#č že bodíky (přenejmenším na vstupu f_modelu) |
|
900 |
|
#č mají konečné souřadnice bez nan a nekonečno. |
|
901 |
|
#č To ale nic neříká o hustotách |
|
902 |
|
mask = np.isfinite(w) |
|
903 |
|
nodes.w = w * mask |
|
904 |
|
#č p_couple není finální pravděpodobnost |
|
905 |
|
#č bodíky ještě musejí projít konvexní obálkou |
|
906 |
|
nodes.p_couple = np.sum(w[mask]) |
|
907 |
|
nodes.score = nodes.p_couple / np.sqrt(np.var(nodes.w)/len(nodes)) |
|
908 |
|
#č a doporučíme (dočasně) k integraci |
|
909 |
|
self._recommend_to_integration(idx) |
|
910 |
|
|
|
911 |
|
|
|
912 |
|
if len(nodes) < nodes.nvar + 1: |
|
913 |
|
# in case of update reset the final score only |
|
914 |
|
nodes.score = -2 |
|
915 |
|
return None |
|
916 |
|
|
|
917 |
|
nodes_model = getattr(nodes, self.model_space) |
|
918 |
|
nodes_pdf = nodes.pdf(self.model_space) |
|
919 |
|
|
|
920 |
|
S_bc = np.cov(nodes_model, rowvar=False, bias=True, aweights=nodes_pdf) |
|
921 |
|
|
|
922 |
|
#č matika |
|
923 |
|
w, v = eig = np.linalg.eig(S_bc) |
|
924 |
|
|
|
925 |
|
if not np.all(w > 0): |
|
926 |
|
nodes.score = -1 |
|
927 |
|
return None |
|
928 |
|
|
|
929 |
|
barycenter = np.average(nodes_model, axis=0, weights=nodes_pdf) |
|
930 |
|
nodes.set_attrs(eig=eig, barycenter=barycenter) |
|
931 |
|
nodes.score = np.trace(S_bc) * np.sqrt(len(nodes)) |
|
932 |
|
|
|
933 |
|
|
|
934 |
|
|
|
935 |
884 |
|
|
936 |
885 |
|
|
937 |
886 |
# nechám s velkým písmenkem |
# nechám s velkým písmenkem |
|
... |
... |
class ContactVoronoi: |
1381 |
1330 |
|
|
1382 |
1331 |
|
|
1383 |
1332 |
## checklist |
## checklist |
1384 |
|
self.couples[tp.couple_indices].append(idx) |
|
|
1333 |
|
self.couples[tp.couple_indices].append(nodes) |
1385 |
1334 |
#self.couples_stats #č nevím |
#self.couples_stats #č nevím |
1386 |
1335 |
#self.mixed_couples #č v této fázi nevíme, |
#self.mixed_couples #č v této fázi nevíme, |
1387 |
1336 |
#self.red_couples #č zda je sada ta nejlepší |
#self.red_couples #č zda je sada ta nejlepší |
|
... |
... |
class ContactVoronoi: |
1493 |
1442 |
|
|
1494 |
1443 |
|
|
1495 |
1444 |
|
|
|
1445 |
|
#č tahle funkce může být volána i skřiňkou |
|
1446 |
|
def estimate_mixed(self, nodes): |
|
1447 |
|
"assigns d1, fv, fw" |
|
1448 |
|
nodes_model = getattr(nodes, self.model_space) |
|
1449 |
|
i, j = nodes.couple |
|
1450 |
|
imask = nodes.imask |
|
1451 |
|
d2 = nodes.d2 |
|
1452 |
|
|
|
1453 |
|
di1 = cdist(nodes_model[imask], np.atleast_2d(self.points[i])).flatten() |
|
1454 |
|
dj1 = cdist(nodes_model[~imask], np.atleast_2d(self.points[j])).flatten() |
|
1455 |
|
|
|
1456 |
|
d1 = np.empty(len(imask)) |
|
1457 |
|
d1[imask] = di1 |
|
1458 |
|
d1[~imask] = dj1 |
|
1459 |
|
#č skříňka bude chtit d1 |
|
1460 |
|
nodes.d1 = d1 |
|
1461 |
|
|
|
1462 |
|
Pdf_i = self.PDF[i] |
|
1463 |
|
Pdf_j = self.PDF[j] |
|
1464 |
|
|
|
1465 |
|
if self.failsi[i]: |
|
1466 |
|
df[imask] = di1 |
|
1467 |
|
df[~imask] = d2[~imask] |
|
1468 |
|
dpdf = df * Pdf_i |
|
1469 |
|
else: |
|
1470 |
|
df[imask] = d2[imask] |
|
1471 |
|
df[~imask] = dj1 |
|
1472 |
|
dpdf = df * Pdf_j |
|
1473 |
|
|
|
1474 |
|
sumpdf = np.empty(len(imask)) |
|
1475 |
|
sumpdf[imask] = di1 * Pdf_i + d2[imask] * Pdf_j |
|
1476 |
|
sumpdf[~imask] = dj1 * Pdf_j + d2[~imask] * Pdf_i |
1496 |
1477 |
|
|
1497 |
|
def estimate_nodes(self, idx): |
|
1498 |
|
pass |
|
|
1478 |
|
w = nodes.w |
|
1479 |
|
|
|
1480 |
|
nodes.fv = w * df / (d1 + d2) |
|
1481 |
|
nodes.fw = w * dpdf / sumpdf |
1499 |
1482 |
|
|
1500 |
1483 |
|
|
1501 |
1484 |
|
|
|
... |
... |
class ContactVoronoi: |
1513 |
1496 |
# set the minimum distance as the standard deviation of IS densisty |
# set the minimum distance as the standard deviation of IS densisty |
1514 |
1497 |
# mindist dělíme dvěma - hranice buňky je na půlcesty |
# mindist dělíme dvěma - hranice buňky je na půlcesty |
1515 |
1498 |
# a dělíme base_r'kem k získání sigmy. Přemejmenším empiricky mně to vychází |
# a dělíme base_r'kem k získání sigmy. Přemejmenším empiricky mně to vychází |
1516 |
|
# (u Norm'u vždy zadaváme směrodatnou odchylku) |
|
1517 |
1499 |
h_i = [stats.norm(sampled_plan_sing[i,j], mindist_sing[i]/2/base_r) for j in range(sample_box.nvar)] |
h_i = [stats.norm(sampled_plan_sing[i,j], mindist_sing[i]/2/base_r) for j in range(sample_box.nvar)] |
1518 |
1500 |
h = f_models.UnCorD(h_i) |
h = f_models.UnCorD(h_i) |
1519 |
1501 |
|
|
|
... |
... |
class ContactVoronoi: |
1753 |
1735 |
|
|
1754 |
1736 |
|
|
1755 |
1737 |
|
|
|
1738 |
|
def get_pf_estimation(self): |
|
1739 |
|
#č zkusme představit, |
|
1740 |
|
#č že někdo bude volat tuhle funkci pořad dokola, |
|
1741 |
|
#č aniž by se přídaly nové vzorky. |
|
1742 |
|
self._update() |
|
1743 |
|
failure, mixed, pfv, pfw = 0 |
|
1744 |
|
|
|
1745 |
|
self.PDF = self.sample_box.pdf(self.model_space) |
|
1746 |
|
|
|
1747 |
|
valid_insides = self._valid_outsides |
|
1748 |
|
is_inside = self.hull.is_inside |
|
1749 |
|
|
|
1750 |
|
for nodes in self.to_integrate.values(): |
|
1751 |
|
if nodes.idx not in valid_insides: |
|
1752 |
|
nodes.inside = mask = is_inside(nodes) |
|
1753 |
|
valid_insides.add(nodes.idx) |
|
1754 |
|
nodes.p_inside = p_inside = np.sum(nodes.w[mask]) |
|
1755 |
|
else: #č kód třídy musí zajistit konzistentnost p_inside, |
|
1756 |
|
try: #č mazajic ji kdyby něco |
|
1757 |
|
p_inside = nodes.p_inside |
|
1758 |
|
except AttributeError: |
|
1759 |
|
nodes.p_inside = p_inside = np.sum(nodes.w[nodes.inside]) |
|
1760 |
|
|
|
1761 |
|
# -1 = 'outside', 0=success, 1=failure, 2=mix |
|
1762 |
|
event_id = nodes.event_id |
|
1763 |
|
if event_id == 1: |
|
1764 |
|
failure += p_inside |
|
1765 |
|
pfv += p_inside |
|
1766 |
|
pfw += p_inside |
|
1767 |
|
|
|
1768 |
|
elif event_id == 2: |
|
1769 |
|
mixed += p_inside |
|
1770 |
|
|
|
1771 |
|
pfv += p_inside |
|
1772 |
|
pfw += p_inside |
|
1773 |
|
else: |
|
1774 |
|
raise ValueError |
|
1775 |
|
|
|
1776 |
|
|
|
1777 |
|
for nodes in self.mixed_couples.values(): |
|
1778 |
|
if nodes.idx not in valid_insides: |
|
1779 |
|
nodes.inside = mask = is_inside(nodes) |
|
1780 |
|
valid_insides.add(nodes.idx) |
|
1781 |
|
nodes.p_inside = p_inside = np.sum(nodes.w[mask]) |
|
1782 |
|
nodes.p_fv = p_fv = np.sum(nodes.fv[mask]) |
|
1783 |
|
nodes.p_fw = p_fw = np.sum(nodes.fw[mask]) |
|
1784 |
|
else: #č kód třídy musí zajistit konzistentnost p_inside, |
|
1785 |
|
try: #č mazajic ji kdyby něco |
|
1786 |
|
p_inside = nodes.p_inside |
|
1787 |
|
p_fv = nodes.p_fv |
|
1788 |
|
p_fw = nodes.p_fw |
|
1789 |
|
except AttributeError: |
|
1790 |
|
nodes.p_inside = p_inside = np.sum(nodes.w[nodes.inside]) |
|
1791 |
|
|
|
1792 |
|
|
|
1793 |
|
|
|
1794 |
|
|
|
1795 |
|
|
|
1796 |
|
|
|
1797 |
|
return failure, mixed, pfv, pfw |
|
1798 |
|
|
|
1799 |
|
# TRI-compatible estimation |
|
1800 |
|
# -1=outside, 0=success, 1=failure, 2=mix |
|
1801 |
|
#č dostaneme vyrovnané odhady Brna-města (-2) a Brna-venkova (-1) |
|
1802 |
|
tri_estimation = dict() |
|
1803 |
|
tri_estimation[-1] = bx.siss.estimations[-1] |
|
1804 |
|
|
|
1805 |
|
# set initial values |
|
1806 |
|
tri_estimation[0] = 0 |
|
1807 |
|
tri_estimation[1] = 0 |
|
1808 |
|
tri_estimation[2] = 0 |
|
1809 |
|
|
|
1810 |
|
#č něco konkretnějšího |
|
1811 |
|
vertex_estimation = 0 |
|
1812 |
|
weighted_vertex_estimation = 0 |
|
1813 |
|
|
|
1814 |
|
pf_inside = bx.siss.estimations[-2] |
|
1815 |
|
if "tri" in dir(bx): |
|
1816 |
|
# let's iterate over all simplices we have |
|
1817 |
|
# (calling code is responsible for the .simplex_stats validity) |
|
1818 |
|
for event_id, simplex_measure, fm, wfm in bx.simplex_stats.values(): |
|
1819 |
|
tri_estimation[event_id] += simplex_measure |
|
1820 |
|
vertex_estimation += fm |
|
1821 |
|
weighted_vertex_estimation += wfm |
|
1822 |
|
|
|
1823 |
|
#č success klidně může být i záporným |
|
1824 |
|
tri_estimation[0] = pf_inside - tri_estimation[1] - tri_estimation[2] |
|
1825 |
|
|
|
1826 |
|
elif np.all(bx.failsi[:len(bx.convex_hull.points)]): |
|
1827 |
|
#č veškerej vnitršek je v poruše |
|
1828 |
|
tri_estimation[1] = pf_inside |
|
1829 |
|
vertex_estimation = pf_inside |
|
1830 |
|
weighted_vertex_estimation = pf_inside |
|
1831 |
|
else: |
|
1832 |
|
#č vnitršek je asi v pořadku |
|
1833 |
|
tri_estimation[0] = pf_inside |
|
1834 |
|
|
|
1835 |
|
return {'TRI_overall_estimations': tri_estimation, \ |
|
1836 |
|
'vertex_estimation' : vertex_estimation, \ |
|
1837 |
|
'weighted_vertex_estimation' : weighted_vertex_estimation} |
|
1838 |
|
|
|
1839 |
|
|
1756 |
1840 |
|
|
1757 |
1841 |
|
|
1758 |
1842 |
|
|