File wellmet/voronoi.py changed (mode: 100644) (index b5cb457..73c7f8d) |
... |
... |
import numpy as np |
5 |
5 |
import numpy.ma as ma |
import numpy.ma as ma |
6 |
6 |
import scipy.stats as stats |
import scipy.stats as stats |
7 |
7 |
|
|
8 |
|
from scipy.spatial import cKDTree |
|
|
8 |
|
from scipy.spatial import KDTree |
9 |
9 |
from scipy.spatial import Delaunay |
from scipy.spatial import Delaunay |
10 |
10 |
from scipy import spatial |
from scipy import spatial |
11 |
11 |
from scipy import interpolate |
from scipy import interpolate |
|
12 |
|
from scipy.optimize import linprog |
12 |
13 |
|
|
13 |
14 |
import collections # for defaultdict |
import collections # for defaultdict |
14 |
15 |
|
|
|
... |
... |
from . import sball |
21 |
22 |
|
|
22 |
23 |
|
|
23 |
24 |
|
|
|
25 |
|
|
|
26 |
|
class ContactSolver: |
|
27 |
|
""" |
|
28 |
|
č Hlavní pointa třídy: |
|
29 |
|
č pokud dva body zvednuté na povrch convexního paraboloidu |
|
30 |
|
č v prostoru ndim+1 (feature space, quadratic kernel) |
|
31 |
|
č lze lineárně separovat (hyperrovinou) od ostatních bodů, |
|
32 |
|
č znamená to, že v původním prostoru |
|
33 |
|
č jejich Voroného buňky budou mít společnou stěnu (kontakt), |
|
34 |
|
č neboli, což je totež, u Delone triangulace by měly společné simplexy |
|
35 |
|
|
|
36 |
|
linprog: |
|
37 |
|
minimize c @ x |
|
38 |
|
such that: |
|
39 |
|
A_ub @ x <= b_ub |
|
40 |
|
A_eq @ x == b_eq |
|
41 |
|
lb <= x <= ub |
|
42 |
|
|
|
43 |
|
č rovnice hyperroviny |
|
44 |
|
H = ax + b = a1x1 + a2x2 + ... + b = 0 |
|
45 |
|
č ačka a bčko jsou pro nás neznamé |
|
46 |
|
č takže máme ndim+1 iksů do tamtoho lineárního solveru |
|
47 |
|
|
|
48 |
|
č nebudeme puntičkařit a pro jednoduchost předepíšeme, |
|
49 |
|
č nechť dva body zájmu leží přímo v hyperrovině |
|
50 |
|
č (bude to hlasít existence rozhraní i v degenerovaných případech |
|
51 |
|
č jako např. tří teček na jedné přímce. Mně to ale nevadí) |
|
52 |
|
č Takže. |
|
53 |
|
|
|
54 |
|
č pro linprog zadáme constrains Ax=b |
|
55 |
|
A = [[ x_1i, x_2i, ..., x_ni, 1], |
|
56 |
|
[ x_1j, x_2j, ..., x_nj, 1]] |
|
57 |
|
b = [0, 0] |
|
58 |
|
|
|
59 |
|
|
|
60 |
|
č Zbytek nechť splnuje nerovnost, tj. |
|
61 |
|
"Convex hull inequalities of the form Ax + b <= 0" |
|
62 |
|
|
|
63 |
|
A = [[ x_1i, x_2i, ..., x_ni, 1], |
|
64 |
|
[ x_1j, x_2j, ..., x_nj, 1]] |
|
65 |
|
b = [0, 0] |
|
66 |
|
|
|
67 |
|
|
|
68 |
|
""" |
|
69 |
|
def __init__(self, points): |
|
70 |
|
nsim, ndim = points.shape |
|
71 |
|
# kind of datascience. feature space, quadratic kernel... |
|
72 |
|
self.lifted_points = np.empty((nsim, ndim + 1)) |
|
73 |
|
self.lifted_points[:, :ndim] = points |
|
74 |
|
self.lifted_points[:, -1] = np.sum(np.square(points), axis=1) |
|
75 |
|
|
|
76 |
|
#č pro linprog zadáme constrains Ax=b |
|
77 |
|
# A = [[ x_1i, x_2i, ..., x_ni, 1], |
|
78 |
|
# [ x_1j, x_2j, ..., x_nj, 1]] |
|
79 |
|
# b = [0, 0] |
|
80 |
|
self.A_eq = np.empty((2, ndim + 2)) |
|
81 |
|
self.A_eq[:, -1] = 1 |
|
82 |
|
self.b_eq = np.zeros(2) |
|
83 |
|
|
|
84 |
|
self.A_ub = np.empty((nsim, ndim + 2)) |
|
85 |
|
self.A_ub[:, :ndim] = points |
|
86 |
|
self.A_ub[:, -1] = np.sum(np.square(points), axis=1) |
|
87 |
|
|
|
88 |
|
def is_couple(self, couple_indices): |
|
89 |
|
i, j = couple_indices |
|
90 |
|
|
|
91 |
|
#č pro linprog zadáme constrains Ax=b |
|
92 |
|
# A = [[ x_1i, x_2i, ..., x_ni, 1], |
|
93 |
|
# [ x_1j, x_2j, ..., x_nj, 1]] |
|
94 |
|
self.A_eq[:, :ndim] = self.lifted_points[[i, j]] |
|
95 |
|
#self.A_eq[:, -1] = 1 |
|
96 |
|
|
|
97 |
|
|
|
98 |
|
|
|
99 |
|
# Ax + b < 0 |
|
100 |
|
|
|
101 |
|
# Ax + b = 0 |
|
102 |
|
|
|
103 |
|
# Ax < b |
|
104 |
|
|
|
105 |
|
c = np.array([-29.0, -45.0, 0.0, 0.0]) |
|
106 |
|
>>> A_ub = np.array([[1.0, -1.0, -3.0, 0.0], |
|
107 |
|
... [-2.0, 3.0, 7.0, -3.0]]) |
|
108 |
|
>>> b_ub = np.array([5.0, -10.0]) |
|
109 |
|
>>> A_eq = np.array([[2.0, 8.0, 1.0, 0.0], |
|
110 |
|
... [4.0, 4.0, 0.0, 1.0]]) |
|
111 |
|
>>> b_eq = np.array([60.0, 60.0]) |
|
112 |
|
>>> x0_bounds = (0, None) |
|
113 |
|
>>> x1_bounds = (0, 5.0) |
|
114 |
|
>>> x2_bounds = (-np.inf, 0.5) # +/- np.inf can be used instead of None |
|
115 |
|
>>> x3_bounds = (-3.0, None) |
|
116 |
|
>>> bounds = [x0_bounds, x1_bounds, x2_bounds, x3_bounds] |
|
117 |
|
>>> result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds) |
|
118 |
|
|
|
119 |
|
|
24 |
120 |
#оӵ кык точкаен Вороной |
#оӵ кык точкаен Вороной |
25 |
121 |
# KickTouchCayenneVoronoi |
# KickTouchCayenneVoronoi |
26 |
122 |
|
|
|
... |
... |
class ContactVoronoi: |
101 |
197 |
nsim = self.sample_box.nsim |
nsim = self.sample_box.nsim |
102 |
198 |
if nsim < 2: |
if nsim < 2: |
103 |
199 |
return None |
return None |
104 |
|
|
|
105 |
|
# # ощень, ощень скромненько |
|
106 |
|
# p_base = 0.3 # chcu, aby behem první iterace p_base tečiček jistě dopadlo do buňky |
|
107 |
|
# nis_base = int((sample_box.nvar+1)/p_base) |
|
108 |
|
# nis = int(max(budget/nsim, nis_base / stats.binom.sf(int(sample_box.nvar+1), nis_base, p_base))) |
|
109 |
200 |
|
|
110 |
201 |
#č já vím, že sample box pokážde failsi přepočítavá |
#č já vím, že sample box pokážde failsi přepočítavá |
111 |
202 |
self.failsi = self.sample_box.failsi |
self.failsi = self.sample_box.failsi |
112 |
203 |
|
|
113 |
204 |
self.PDF = self.sample_box.pdf(self.model_space) |
self.PDF = self.sample_box.pdf(self.model_space) |
114 |
205 |
|
|
115 |
|
#č zde provadím rozdělení na prostor, ve kterém vzorkujem |
|
116 |
|
#č a prostor "modelu", vô ktôrom, v podstatě, měříme vzdaleností |
|
117 |
|
sampled_plan_model = getattr(self.sample_box, model_space) |
|
118 |
|
self.tree = cKDTree(sampled_plan_model) |
|
119 |
|
|
|
|
206 |
|
self.sampled_plan_model = getattr(self.sample_box, self.model_space) |
|
207 |
|
self.tree = KDTree(sampled_plan_model) |
120 |
208 |
|
|
121 |
|
if self.sampling_space is None: |
|
122 |
|
sampling_space = model_space |
|
123 |
|
# sing like sampling |
|
124 |
|
sampled_plan_sing = sampled_plan_model |
|
125 |
|
self.tree_sampling = self.tree |
|
126 |
|
else: |
|
127 |
|
sampled_plan_sing = getattr(self.sample_box, self.sampling_space) |
|
128 |
|
# narychlo, moc to nepomůže, neboť asi po 16 vzorcích počítá brut forsem |
|
129 |
|
self.tree_sampling = cKDTree(sampled_plan_sing, compact_nodes=False, balanced_tree=False) |
|
130 |
|
|
|
131 |
|
# find distance to the nearest sampling point (from all points) |
|
132 |
|
#č pro sampling bych vždy použival p_norm=2 |
|
133 |
|
dd2, ii2 = tree_sampling.query(sampled_plan_sing, k=[2], p=2) |
|
134 |
|
mindist_sing = dd2.flatten() |
|
135 |
209 |
|
|
136 |
210 |
|
|
137 |
211 |
#č 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: |
194 |
268 |
def onboard_couple(self, couple_indices): |
def onboard_couple(self, couple_indices): |
195 |
269 |
#č (i, j)? poprvé ta čísla vidíme |
#č (i, j)? poprvé ta čísla vidíme |
196 |
270 |
i, j = couple_indices |
i, j = couple_indices |
197 |
|
if i in self.failure_points: |
|
198 |
|
if j in self.failure_points: |
|
199 |
|
#č červený kontakt |
|
200 |
|
self.integrate_couple(i) |
|
201 |
|
else: |
|
202 |
|
#č žlutý kontakt |
|
203 |
|
self.integrate_couple(i) |
|
204 |
|
else: |
|
205 |
|
|
|
|
271 |
|
|
|
272 |
|
contact_node_model = np.sum(self.sampled_plan_model[[i, j]], axis=0) |
|
273 |
|
np.divide(contact_node_model, 2, out=contact_node_model) |
|
274 |
|
|
|
275 |
|
|
206 |
276 |
self.direct_contacts[(i, j)] = None |
self.direct_contacts[(i, j)] = None |
207 |
277 |
self.nodes[(i, j)] = None |
self.nodes[(i, j)] = None |
208 |
278 |
|
|
209 |
279 |
|
|
|
280 |
|
dd, ii = tree.query(h_plan_model_part, k=1, p=p_norm) |
|
281 |
|
|
|
282 |
|
# find distance to the nearest sampling point (from all points) |
|
283 |
|
dd2, ii2 = tree_sampling.query(sampled_plan_sing, k=[2], p=2) |
|
284 |
|
mindist_sing = dd2.flatten() |
|
285 |
|
|
210 |
286 |
|
|
211 |
287 |
|
|
212 |
288 |
|
|
|
... |
... |
class ContactVoronoi: |
220 |
296 |
self.onboard_couple((i, j)) |
self.onboard_couple((i, j)) |
221 |
297 |
|
|
222 |
298 |
|
|
|
299 |
|
#č i-té indexy jsou ty čerstvé, přes ně iterujeme |
|
300 |
|
#č j-té - ty bežné, protí ním rozhodujeme |
|
301 |
|
if j < self._nsim: |
|
302 |
|
self._indices_to_update.add(j) |
|
303 |
|
|
|
304 |
|
|
223 |
305 |
|
|
224 |
306 |
def integrate_couple(self, couple_indices): |
def integrate_couple(self, couple_indices): |
225 |
307 |
"""č máme nový manželský páreček |
"""č máme nový manželský páreček |
|
... |
... |
class ContactVoronoi: |
246 |
328 |
h_plan_model = getattr(h_plan, model_space) |
h_plan_model = getattr(h_plan, model_space) |
247 |
329 |
dd, ii = tree.query(h_plan_model, k=1, p=p_norm) |
dd, ii = tree.query(h_plan_model, k=1, p=p_norm) |
248 |
330 |
|
|
|
331 |
|
|
249 |
332 |
# nechám s velkým písmenkem |
# nechám s velkým písmenkem |
250 |
333 |
Vor_mask = ii==i |
Vor_mask = ii==i |
251 |
334 |
h_plan_model_ma = h_plan_model[Vor_mask] |
h_plan_model_ma = h_plan_model[Vor_mask] |