File wellmet/voronoi.py changed (mode: 100644) (index 73c7f8d..7a39675) |
... |
... |
class ContactSolver: |
57 |
57 |
b = [0, 0] |
b = [0, 0] |
58 |
58 |
|
|
59 |
59 |
|
|
60 |
|
č Zbytek nechť splnuje nerovnost, tj. |
|
|
60 |
|
č Zbytek musí splňovat nerovnost, tj. |
61 |
61 |
"Convex hull inequalities of the form Ax + b <= 0" |
"Convex hull inequalities of the form Ax + b <= 0" |
62 |
62 |
|
|
63 |
|
A = [[ x_1i, x_2i, ..., x_ni, 1], |
|
64 |
|
[ x_1j, x_2j, ..., x_nj, 1]] |
|
65 |
|
b = [0, 0] |
|
66 |
|
|
|
|
63 |
|
č Neboli Ax <= b v termínech linprog |
|
64 |
|
č díky "menší nebo rovno" |
|
65 |
|
č nemusíme plnou matici ani maskovat |
67 |
66 |
|
|
|
67 |
|
č ačka hyperroviny zadavají jednotkový vektor normály, |
|
68 |
|
č takže leží v mezích -1 < a < 1 |
|
69 |
|
č bčko lze omezit poloosou |
68 |
70 |
""" |
""" |
69 |
71 |
def __init__(self, points): |
def __init__(self, points): |
70 |
72 |
nsim, ndim = points.shape |
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 |
73 |
|
|
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) |
|
|
74 |
|
self.A = A = np.empty((nsim, ndim + 2)) |
|
75 |
|
A[:, :ndim] = points |
|
76 |
|
# kind of datascience. feature space, quadratic kernel... |
|
77 |
|
#č dáme to trochu niž (ten "ndim"), abychom se vyhli triviálnímu řešení |
|
78 |
|
A[:, -2] = np.sum(np.square(points)-ndim, axis=1) |
|
79 |
|
A[:, -1] = 1 |
|
80 |
|
|
|
81 |
|
#č žšmaria, alokovat nuly.. |
|
82 |
|
self.b_ub = np.zeros(nsim) |
|
83 |
|
self.c = np.zeros(ndim + 2) |
|
84 |
|
self.b_eq = self.c[:2] |
|
85 |
|
self.bounds = [(-1, 1) for __ in range(ndim + 1)] |
|
86 |
|
#č skoro-nulové b-čko může být legální |
|
87 |
|
#č ale musíme vyhnout triviálnímu řešení |
|
88 |
|
#č nevím co s tím |
|
89 |
|
self.bounds.append((None, -0.3)) |
87 |
90 |
|
|
88 |
91 |
def is_couple(self, couple_indices): |
def is_couple(self, couple_indices): |
89 |
92 |
i, j = couple_indices |
i, j = couple_indices |
|
... |
... |
class ContactSolver: |
91 |
94 |
#č pro linprog zadáme constrains Ax=b |
#č pro linprog zadáme constrains Ax=b |
92 |
95 |
# A = [[ x_1i, x_2i, ..., x_ni, 1], |
# A = [[ x_1i, x_2i, ..., x_ni, 1], |
93 |
96 |
# [ x_1j, x_2j, ..., x_nj, 1]] |
# [ 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 |
|
|
97 |
|
A_eq = self.A[[i, j]] |
100 |
98 |
|
|
101 |
|
# Ax + b = 0 |
|
|
99 |
|
#č nevím, zda omezení (bounds) zrychlí, nebo zpomalí běh řešiče |
|
100 |
|
result = linprog(self.c, A_ub=self.A, b_ub=self.b_ub, A_eq=A_eq,\ |
|
101 |
|
b_eq=self.b_eq, bounds=self.bounds) |
102 |
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) |
|
|
103 |
|
return result.success |
|
104 |
|
|
118 |
105 |
|
|
119 |
106 |
|
|
120 |
107 |
#оӵ кык точкаен Вороной |
#оӵ кык точкаен Вороной |
|
... |
... |
class ContactVoronoi: |
154 |
141 |
č 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 |
155 |
142 |
""" |
""" |
156 |
143 |
def __init__(self, sample_box, hull, model_space='G', \ |
def __init__(self, sample_box, hull, model_space='G', \ |
157 |
|
p_norm=2, p_base=0.5, ns=1000, auto_update=True, |
|
158 |
|
sample_callback=None) |
|
|
144 |
|
p_norm=2, p_base=0.5, ns=1000, auto_update=True, \ |
|
145 |
|
sample_callback=None): |
159 |
146 |
self.sample_box = sample_box |
self.sample_box = sample_box |
160 |
147 |
self.hull = hull |
self.hull = hull |
161 |
148 |
self.model_space = model_space |
self.model_space = model_space |
|
... |
... |
class ContactVoronoi: |
247 |
234 |
#č u prvního indexu musíme dodržovat pořadí procházky |
#č u prvního indexu musíme dodržovat pořadí procházky |
248 |
235 |
for i in range(self._nsim - 1, 0, -1): |
for i in range(self._nsim - 1, 0, -1): |
249 |
236 |
if i in self._indices_to_update: |
if i in self._indices_to_update: |
250 |
|
self._indices_to_update.remove(i)� |
|
|
237 |
|
self._indices_to_update.remove(i) |
251 |
238 |
for j in self._indices_to_update: |
for j in self._indices_to_update: |
252 |
239 |
if ((i, j) in self.mixed_couples) or ((i, j) in self.red_couples): |
if ((i, j) in self.mixed_couples) or ((i, j) in self.red_couples): |
253 |
240 |
self.update_couple((i, j)) |
self.update_couple((i, j)) |
|
... |
... |
class ContactVoronoi: |
273 |
260 |
np.divide(contact_node_model, 2, out=contact_node_model) |
np.divide(contact_node_model, 2, out=contact_node_model) |
274 |
261 |
|
|
275 |
262 |
|
|
276 |
|
self.direct_contacts[(i, j)] = None |
|
277 |
|
self.nodes[(i, j)] = None |
|
|
263 |
|
# self.direct_contacts[(i, j)] = None |
|
264 |
|
# self.nodes[(i, j)] = None |
278 |
265 |
|
|
279 |
266 |
|
|
280 |
267 |
dd, ii = tree.query(h_plan_model_part, k=1, p=p_norm) |
dd, ii = tree.query(h_plan_model_part, k=1, p=p_norm) |
|
... |
... |
class ContactVoronoi: |
287 |
274 |
|
|
288 |
275 |
|
|
289 |
276 |
#č není znamená, že nikdy ani nebylo |
#č není znamená, že nikdy ani nebylo |
290 |
|
if ((i, j) in self.nodes) or ((i, j) in self.direct_contacts): |
|
291 |
|
#č None znamená, že kontakt je pro nás pase |
|
292 |
|
if self.nodes[(i, j)] is not None: |
|
293 |
|
self.update_couple((i, j)) |
|
294 |
|
|
|
295 |
|
else: #č (i, j)? poprvé ta čísla vidíme |
|
296 |
|
self.onboard_couple((i, j)) |
|
|
277 |
|
# if ((i, j) in self.nodes) or ((i, j) in self.direct_contacts): |
|
278 |
|
# #č None znamená, že kontakt je pro nás pase |
|
279 |
|
# if self.nodes[(i, j)] is not None: |
|
280 |
|
# self.update_couple((i, j)) |
|
281 |
|
# |
|
282 |
|
# else: #č (i, j)? poprvé ta čísla vidíme |
|
283 |
|
# self.onboard_couple((i, j)) |
297 |
284 |
|
|
298 |
285 |
|
|
299 |
286 |
#č i-té indexy jsou ty čerstvé, přes ně iterujeme |
#č i-té indexy jsou ty čerstvé, přes ně iterujeme |