File wellmet/voronoi.py changed (mode: 100644) (index 5f798fd..8624ad2) |
... |
... |
class ConvexSolver: |
111 |
111 |
return is_line_convex(self.lifted_points, couple_indices) |
return is_line_convex(self.lifted_points, couple_indices) |
112 |
112 |
|
|
113 |
113 |
|
|
|
114 |
|
#č Odbočka z odbočky |
|
115 |
|
class ConvexSpline: |
|
116 |
|
""" |
|
117 |
|
č Hlavní pointa CS tříd: |
|
118 |
|
č pokud dva body zvednuté na povrch convexního paraboloidu |
|
119 |
|
č v prostoru ndim+1 (feature space, quadratic kernel) |
|
120 |
|
č lze lineárně separovat (hyperrovinou) od ostatních bodů, |
|
121 |
|
č znamená to, že v původním prostoru |
|
122 |
|
č jejich Voroného buňky budou mít společnou stěnu (kontakt), |
|
123 |
|
č neboli, což je totež, u Delone triangulace by měly společné simplexy |
|
124 |
|
|
|
125 |
|
č Hlavní pointa ConvexSpline třídy: |
|
126 |
|
č Využit navíc geometrických informací, které už předem známé: |
|
127 |
|
č 1. Známe souřadnice vzorků. |
|
128 |
|
č 2. Víme, že přímka mezí těmi dvěma vzorky leží v hyperrovině |
|
129 |
|
č 3. Vždyť to my zvedáme na povrch convexního paraboloidu! |
|
130 |
|
č Můžeme tedy v každém její bodě najit tečnou hyperrovinu. |
|
131 |
|
|
|
132 |
|
|
|
133 |
|
č nebudeme puntičkařit a pro jednoduchost předepíšeme, |
|
134 |
|
č nechť dva body zájmu leží přímo v hyperrovině |
|
135 |
|
č (bude to hlasít existence rozhraní i v degenerovaných případech |
|
136 |
|
č jako např. tří teček na jedné přímce. Mně to ale nevadí) |
|
137 |
|
""" |
|
138 |
|
def __init__(self, points): |
|
139 |
|
nsim, ndim = points.shape |
|
140 |
|
|
|
141 |
|
self.lifted_points = np.empty((nsim, ndim + 1)) |
|
142 |
|
self.lifted_points[:, :ndim] = points |
|
143 |
|
# kind of datascience. feature space, quadratic kernel... |
|
144 |
|
self.lifted_points[:, -1] = np.sum(np.square(points), axis=1) |
|
145 |
|
|
|
146 |
|
|
|
147 |
|
def is_couple(self, couple_indices): |
|
148 |
|
i, j = couple_indices |
|
149 |
|
|
|
150 |
|
X = self.lifted_points |
|
151 |
|
__nsim, ndim = X.shape |
|
152 |
|
|
|
153 |
|
|
|
154 |
|
|
|
155 |
|
|
|
156 |
|
basis = np.random.random((ndim, ndim)) |
|
157 |
|
#č první vektor musí být zadán přímkou mezí vzorky |
|
158 |
|
#č jinak se posype náš předpoklad, že leží v hyperrovině |
|
159 |
|
#č a žádné výsledky "za", "před" hyperrovinou nebudou nic znamenat |
|
160 |
|
basis[:, 0] = X[j] - X[i] #č QR rozklad jede po sloupcich |
|
161 |
|
|
|
162 |
|
#č jako indice zkusme použit bod na usečce uprostřed mezí vzorky |
|
163 |
|
half_point = np.sum(self.lifted_points[[i, j]], axis=0) |
|
164 |
|
np.divide(half_point, 2, out=half_point) |
|
165 |
|
|
|
166 |
|
#č všechy souřadnice jsou dány radius-vektorem od středu, |
|
167 |
|
#č ale poslední souřadnice je to naše zvednutí, |
|
168 |
|
#č kde bychom mohli zkusit dát korrektnější směr. |
|
169 |
|
#č Mně humpolackými uvahami o tečně paraboly |
|
170 |
|
#č na caru paríru vyšlo něco jako |
|
171 |
|
#č že stačí dát poslední složku 0,5. |
|
172 |
|
#č Jakože čím je roloměr je větší, |
|
173 |
|
#č tím je "svislá" složka automaticky měnší. |
|
174 |
|
#č Jakože netřeba ani normalizovat, ani nic "složitě" počítat |
|
175 |
|
half_point[-1] = -0.5 |
|
176 |
|
|
|
177 |
|
#č ten náš radius-vektor není ortogonální |
|
178 |
|
#č k přímce, na které leží ty naši dva vzorky |
|
179 |
|
#č QR rozklad je ortogonalizuje |
|
180 |
|
#č a vygeneruje další ortogonalní vektory |
|
181 |
|
basis[:, 1] = half_point |
|
182 |
|
basis, __ = np.linalg.qr(basis) |
|
183 |
|
|
|
184 |
|
|
|
185 |
|
#č vytahneme náš ortogonalizovaný odhad |
|
186 |
|
#č normálního vektoru |
|
187 |
|
a = basis[:, 1] |
|
188 |
|
##č nejak tuším, že v poslední dimenzi |
|
189 |
|
##č znaménko normál musí být jednotné |
|
190 |
|
a = a * -np.sign(a[-1]) |
|
191 |
|
b = X @ a |
|
192 |
|
|
|
193 |
|
if np.max(b) - np.max(b[[i, j]]) < 1e-7: |
|
194 |
|
#č bylo to myšleno jako jakési zoufalství |
|
195 |
|
#č ale funguje překvapivě dobře |
|
196 |
|
return True |
|
197 |
|
else: |
|
198 |
|
idx = np.argmax(b) |
|
199 |
|
#č nahradíme v jedničce naš odhad normálního vektoru |
|
200 |
|
#č nalezenou překažkou. |
|
201 |
|
#č Zbytek bazí je jíž byl ortogonální |
|
202 |
|
#č k té naši pomyšlené normále |
|
203 |
|
basis[:, 1] = X[idx] - X[i] |
|
204 |
|
|
|
205 |
|
#č Ve zbytku pokračujeme jako vždycky. |
|
206 |
|
#č Pořad ndim pokusu u False párečku |
|
207 |
|
for __ in range(ndim): |
|
208 |
|
basis, __ = np.linalg.qr(basis) |
|
209 |
|
|
|
210 |
|
# get constrain |
|
211 |
|
a = basis[:, -1] |
|
212 |
|
#č nejak tuším, že v poslední dimenzi |
|
213 |
|
#č znaménko normál musí být jednotné |
|
214 |
|
a = a * -np.sign(a[-1]) |
|
215 |
|
b = X @ a |
|
216 |
|
if np.max(b) - np.max(b[[i, j]]) < 1e-7: |
|
217 |
|
return True |
|
218 |
|
else: |
|
219 |
|
idx = np.argmax(b) |
|
220 |
|
basis[:, 2:] = basis[:, 1:-1] |
|
221 |
|
basis[:, 1] = X[idx] - X[i] |
|
222 |
|
|
|
223 |
|
basis, __ = np.linalg.qr(basis) |
|
224 |
|
a = basis[:, -1] |
|
225 |
|
#č nejak tuším, že v poslední dimenzi |
|
226 |
|
#č znaménko normál musí být jednotné |
|
227 |
|
a = a * -np.sign(a[-1]) |
|
228 |
|
b = X @ a |
|
229 |
|
return np.max(b) - np.max(b[[i, j]]) < 1e-7 |
114 |
230 |
|
|
115 |
231 |
|
|
116 |
232 |
|
|
|
... |
... |
class ContactVoronoi: |
262 |
378 |
č Neplatí: |
č Neplatí: |
263 |
379 |
č 2. (Při přidani vzorku "c"). Jakože neexistují-li kontakty mezi vzorky "c" a "a", |
č 2. (Při přidani vzorku "c"). Jakože neexistují-li kontakty mezi vzorky "c" a "a", |
264 |
380 |
č a mezi "c" a "b", tak přídaní "c" nemůže kontakt mezi "a" a "b" ovlivnit. |
č a mezi "c" a "b", tak přídaní "c" nemůže kontakt mezi "a" a "b" ovlivnit. |
|
381 |
|
|
|
382 |
|
|
|
383 |
|
č Poznámky k adaptivnímu vzorkování: |
|
384 |
|
č Pro vzorkování alike potřebujeme CoV, |
|
385 |
|
č kde važení hustotámi je relativné, |
|
386 |
|
č násobení dvakrát aweights nezmění výslednou kovarianční matici. |
|
387 |
|
č Pro rozhodování, která sada nejlepé zachytila "buňku" - |
|
388 |
|
č něco jako matici setrvačnosti, u které by násobění "hmotností" |
|
389 |
|
č konstantou zvětšovalo by celkovou setrvačnost systému. |
|
390 |
|
č Tušíme, že (obecně) dobře zachycená "buňka" by měla výkazovat |
|
391 |
|
č největší setrvačnost, ale i taky prostě hmotnost (pravděpodobnostní obsah). |
|
392 |
|
č Mohli bychom zůstat u CoV pro účely adaptivního vzorkování |
|
393 |
|
č a volít nejlepší sadu podle podílu pravděpodobnostního obsahu |
|
394 |
|
č k rozptylu estimatoru |
|
395 |
|
|
|
396 |
|
|
|
397 |
|
|
|
398 |
|
č co máme zadavat jako aweight u np.cov? |
|
399 |
|
č pdf bodíku? Nebo jejich vahy? |
|
400 |
|
|
|
401 |
|
č jak počítat skore? čím máme dělit pravděpodobnostní obsah? |
|
402 |
|
č 1. sum((w_i - w_mean)**2) / n To je rozptyl w |
|
403 |
|
č 2. sum((w_i - w_mean)**2) / n**2 To je rozptyl odhadu w_mean |
|
404 |
|
č 3. sqrt(sum((w_i - w_mean)**2)) / n Směrodatná odhylka odhadu w_mean |
265 |
405 |
""" |
""" |
266 |
406 |
def __init__(self, sample_box, hull, model_space='G', \ |
def __init__(self, sample_box, hull, model_space='G', \ |
267 |
407 |
p_norm=2, ns=1000, p_base=0.5, auto_update=True, \ |
p_norm=2, ns=1000, p_base=0.5, auto_update=True, \ |
|
... |
... |
class ContactVoronoi: |
359 |
499 |
#č já vím, že sample box pokážde failsi přepočítavá |
#č já vím, že sample box pokážde failsi přepočítavá |
360 |
500 |
self.failsi = failsi = self.sample_box.failsi |
self.failsi = failsi = self.sample_box.failsi |
361 |
501 |
|
|
362 |
|
self.CS = ConvexSolver(sampled_plan_model) |
|
|
502 |
|
self.CS = ConvexSpline(sampled_plan_model) |
363 |
503 |
self.PDF = self.sample_box.pdf(self.model_space) |
self.PDF = self.sample_box.pdf(self.model_space) |
364 |
504 |
|
|
365 |
505 |
assert len(self._indices_to_update) == 0 |
assert len(self._indices_to_update) == 0 |
|
... |
... |
class ContactVoronoi: |
443 |
583 |
# #č a na ně vykašle |
# #č a na ně vykašle |
444 |
584 |
# self._indices_to_update.add((max(j,k), min(j,k))) |
# self._indices_to_update.add((max(j,k), min(j,k))) |
445 |
585 |
|
|
446 |
|
def _recommend_to_integrate(self, nodes_idx): |
|
|
586 |
|
def _recommend_to_integration(self, nodes_idx): |
447 |
587 |
#č Má cenu integrovat pouze nejreprezentativnější sadu bodů. |
#č Má cenu integrovat pouze nejreprezentativnější sadu bodů. |
448 |
588 |
#č Špatné sady je velmi obtižné správně započítavat |
#č Špatné sady je velmi obtižné správně započítavat |
449 |
589 |
#č aby nezhoršovaly výsledek. |
#č aby nezhoršovaly výsledek. |
|
... |
... |
class ContactVoronoi: |
481 |
621 |
np.divide(contact_node_model, 2, out=contact_node_model) |
np.divide(contact_node_model, 2, out=contact_node_model) |
482 |
622 |
|
|
483 |
623 |
#č bodík nebude mít vahu - nechť si s tím integrovačka poradí! |
#č bodík nebude mít vahu - nechť si s tím integrovačka poradí! |
|
624 |
|
#č (ten bodík by neměl být omylem doporučen k integraci) |
484 |
625 |
idx, dd, ii, mask = self._store_masked(contact_node_model, couple_indices, event_id) |
idx, dd, ii, mask = self._store_masked(contact_node_model, couple_indices, event_id) |
485 |
626 |
|
|
486 |
627 |
#č kontakt neexistuje, netřeba nic robiť, j-ko i-čkem není dotčen |
#č kontakt neexistuje, netřeba nic robiť, j-ko i-čkem není dotčen |
|
... |
... |
class ContactVoronoi: |
516 |
657 |
#č uspěch. Dodáme bodíkům vahy |
#č uspěch. Dodáme bodíkům vahy |
517 |
658 |
nodes = self.nodes[idx] |
nodes = self.nodes[idx] |
518 |
659 |
w = nodes.pdf(self.model_space) / pdf[mask] / self.ns |
w = nodes.pdf(self.model_space) / pdf[mask] / self.ns |
519 |
|
nodes.w = w |
|
|
660 |
|
#č Tamta filtrace stromem zaručuje, |
|
661 |
|
#č že bodíky (přenejmenším na vstupu f_modelu) |
|
662 |
|
#č mají konečné souřadnice bez nan a nekonečno. |
|
663 |
|
#č To ale nic neříká o hustotách |
|
664 |
|
mask = np.isfinite(w) |
|
665 |
|
nodes.w = w * mask |
|
666 |
|
nodes.p_couple = np.sum(w[mask]) |
520 |
667 |
#č a doporučíme (dočasně) k integraci |
#č a doporučíme (dočasně) k integraci |
521 |
|
self._recommend_to_integrate(idx) |
|
|
668 |
|
self._recommend_to_integration(idx) |
522 |
669 |
|
|
523 |
670 |
if len(nodes) > ndim: |
if len(nodes) > ndim: |
524 |
671 |
#č a spustíme adaptivní IS vzorkování |
#č a spustíme adaptivní IS vzorkování |
|
... |
... |
class ContactVoronoi: |
781 |
928 |
|
|
782 |
929 |
|
|
783 |
930 |
def score(self, nodes): |
def score(self, nodes): |
784 |
|
|
|
|
931 |
|
""" |
|
932 |
|
č z |
|
933 |
|
""" |
785 |
934 |
if len(nodes) < nodes.nvar + 1: |
if len(nodes) < nodes.nvar + 1: |
786 |
935 |
# in case of update reset the final score only |
# in case of update reset the final score only |
787 |
936 |
nodes.score = -2 |
nodes.score = -2 |
|
... |
... |
class ContactVoronoi: |
789 |
938 |
|
|
790 |
939 |
nodes_model = getattr(nodes, self.model_space) |
nodes_model = getattr(nodes, self.model_space) |
791 |
940 |
nodes_pdf = nodes.pdf(self.model_space) |
nodes_pdf = nodes.pdf(self.model_space) |
|
941 |
|
|
792 |
942 |
S_bc = np.cov(nodes_model, rowvar=False, bias=True, aweights=nodes_pdf) |
S_bc = np.cov(nodes_model, rowvar=False, bias=True, aweights=nodes_pdf) |
793 |
943 |
|
|
794 |
944 |
#č matika |
#č matika |
|
... |
... |
class ContactVoronoi: |
805 |
955 |
|
|
806 |
956 |
|
|
807 |
957 |
|
|
808 |
|
def sample_couple(self, couple_indices, event_id, center, r): |
|
809 |
|
"""č máme nový manželský páreček |
|
810 |
|
оӵ Выль кузо. Мар меда кароно? |
|
|
958 |
|
def sample_alike(self, nodes): |
|
959 |
|
""" |
811 |
960 |
č zde jenom vzorkujeme; |
č zde jenom vzorkujeme; |
812 |
961 |
č integrování, vyhodnocování je jinde a později! |
č integrování, vyhodnocování je jinde a později! |
813 |
962 |
""" |
""" |
|
... |
... |
class ContactVoronoi: |
1094 |
1243 |
f = self.sample_box.f_model |
f = self.sample_box.f_model |
1095 |
1244 |
#č tady je to super, ťažké datové struktury f-modelu |
#č tady je to super, ťažké datové struktury f-modelu |
1096 |
1245 |
#č vytvaříme na jíž vyfiltrovaných datéch |
#č vytvaříme na jíž vyfiltrovaných datéch |
|
1246 |
|
#č Tahle filtrace stromem zaručuje, |
|
1247 |
|
#č že bodíky (přenejmenším na vstupu f_modelu) |
|
1248 |
|
#č mají konečné souřadnice bez nan a nekonečno. |
|
1249 |
|
#č To ale nic neříká o hustotách |
1097 |
1250 |
nodes = f.new_sample(nodes_model[mask], self.model_space) |
nodes = f.new_sample(nodes_model[mask], self.model_space) |
1098 |
1251 |
couple_stats = {'event_id':event_id, 'couple':couple_indices} |
couple_stats = {'event_id':event_id, 'couple':couple_indices} |
1099 |
1252 |
nodes = CandyNodes(nodes, couple_stats) |
nodes = CandyNodes(nodes, couple_stats) |
1100 |
1253 |
idx = self.nodes.add(nodes) |
idx = self.nodes.add(nodes) |
|
1254 |
|
nodes.idx = idx |
1101 |
1255 |
|
|
1102 |
1256 |
#č červené kontakty prostě rovnou integrujeme |
#č červené kontakty prostě rovnou integrujeme |
1103 |
1257 |
#č takže netřeba pro ně ukladat zbytečná data |
#č takže netřeba pro ně ukladat zbytečná data |