File wellmet/wireframe.py changed (mode: 100644) (index 35bf320..4a80e37) |
3 |
3 |
|
|
4 |
4 |
import numpy as np |
import numpy as np |
5 |
5 |
from scipy.spatial import Delaunay, KDTree |
from scipy.spatial import Delaunay, KDTree |
|
6 |
|
from scipy.optimize import linprog |
6 |
7 |
|
|
7 |
8 |
|
|
8 |
9 |
|
|
|
... |
... |
class LocalizedHull: |
729 |
730 |
|
|
730 |
731 |
|
|
731 |
732 |
|
|
732 |
|
#class ContactSolver: |
|
733 |
|
# """ |
|
734 |
|
# č Hlavní pointa třídy: |
|
735 |
|
# č pokud dva body zvednuté na povrch convexního paraboloidu |
|
736 |
|
# č v prostoru ndim+1 (feature space, quadratic kernel) |
|
737 |
|
# č lze lineárně separovat (hyperrovinou) od ostatních bodů, |
|
738 |
|
# č znamená to, že v původním prostoru |
|
739 |
|
# č jejich Voroného buňky budou mít společnou stěnu (kontakt), |
|
740 |
|
# č neboli, což je totež, u Delone triangulace by měly společné simplexy |
|
741 |
|
# |
|
742 |
|
# linprog: |
|
743 |
|
# minimize c @ x |
|
744 |
|
# such that: |
|
745 |
|
# A_ub @ x <= b_ub |
|
746 |
|
# A_eq @ x == b_eq |
|
747 |
|
# lb <= x <= ub |
|
748 |
|
# |
|
749 |
|
# č rovnice hyperroviny |
|
750 |
|
# H = ax + b = a1x1 + a2x2 + ... + b = 0 |
|
751 |
|
# č ačka a bčko jsou pro nás neznamé |
|
752 |
|
# č takže máme ndim+1 iksů do tamtoho lineárního solveru |
|
753 |
|
# |
|
754 |
|
# č nebudeme puntičkařit a pro jednoduchost předepíšeme, |
|
755 |
|
# č nechť dva body zájmu leží přímo v hyperrovině |
|
756 |
|
# č (bude to hlasít existence rozhraní i v degenerovaných případech |
|
757 |
|
# č jako např. tří teček na jedné přímce. Mně to ale nevadí) |
|
758 |
|
# č Takže. |
|
759 |
|
# |
|
760 |
|
# č pro linprog zadáme constrains Ax=b |
|
761 |
|
# A = [[ x_1i, x_2i, ..., x_ni, 1], |
|
762 |
|
# [ x_1j, x_2j, ..., x_nj, 1]] |
|
763 |
|
# b = [0, 0] |
|
764 |
|
# |
|
765 |
|
# |
|
766 |
|
# č Zbytek musí splňovat nerovnost, tj. |
|
767 |
|
# "Convex hull inequalities of the form Ax + b <= 0" |
|
768 |
|
# |
|
769 |
|
# č Neboli Ax <= b v termínech linprog |
|
770 |
|
# č díky "menší nebo rovno" |
|
771 |
|
# č nemusíme plnou matici ani maskovat |
|
772 |
|
# |
|
773 |
|
# č ačka hyperroviny zadavají jednotkový vektor normály, |
|
774 |
|
# č takže leží v mezích -1 < a < 1 |
|
775 |
|
# č bčko lze omezit poloosou |
|
776 |
|
# """ |
|
777 |
|
# def __init__(self, points): |
|
778 |
|
# nsim, ndim = points.shape |
|
779 |
|
# |
|
780 |
|
# self.A = A = np.empty((nsim, ndim + 2)) |
|
781 |
|
# A[:, :ndim] = points |
|
782 |
|
# # kind of datascience. feature space, quadratic kernel... |
|
783 |
|
# #č dáme to trochu niž (ten "ndim"), abychom se vyhli triviálnímu řešení |
|
784 |
|
# A[:, -2] = np.sum(np.square(points), axis=1)-1000 |
|
785 |
|
# A[:, -1] = 1 |
|
786 |
|
# |
|
787 |
|
# self.lifted_points = A[:, :ndim+1] |
|
788 |
|
# |
|
789 |
|
# #č žšmaria, alokovat nuly.. |
|
790 |
|
# self.b_ub = np.atleast_1d(np.zeros(ndim-2)) |
|
791 |
|
# self.c = np.zeros(ndim + 2) |
|
792 |
|
# self.b_eq = self.c[:2] |
|
793 |
|
# self.bounds = [(-1, 1) for __ in range(ndim + 1)] |
|
794 |
|
# #č skoro-nulové b-čko může být legální |
|
795 |
|
# #č ale musíme vyhnout triviálnímu řešení |
|
796 |
|
# #č nevím co s tím |
|
797 |
|
# self.bounds.append((None, -0.2)) |
|
798 |
|
# |
|
799 |
|
# def is_couple(self, couple_indices): |
|
800 |
|
# i, j = couple_indices |
|
801 |
|
# |
|
802 |
|
# #č pro linprog zadáme constrains Ax=b |
|
803 |
|
# # A = [[ x_1i, x_2i, ..., x_ni, 1], |
|
804 |
|
# # [ x_1j, x_2j, ..., x_nj, 1]] |
|
805 |
|
# A_eq = self.A[[i, j]] |
|
806 |
|
# |
|
807 |
|
# mask = np.ones(len(self.A), dtype=np.bool8) |
|
808 |
|
# mask[[i, j]] = False |
|
809 |
|
# A_ub = np.atleast_2d(self.A[mask]) |
|
810 |
|
# |
|
811 |
|
# result = linprog(self.c, A_ub=A_ub, b_ub=self.b_ub, A_eq=A_eq,\ |
|
812 |
|
# b_eq=self.b_eq, bounds=self.bounds) |
|
813 |
|
# |
|
814 |
|
# return result.success |
|
|
733 |
|
|
|
734 |
|
# |
|
735 |
|
# Voronoi adjacency by linear programming |
|
736 |
|
# |
|
737 |
|
|
|
738 |
|
|
|
739 |
|
|
|
740 |
|
class LinearSolver: |
|
741 |
|
""" |
|
742 |
|
linprog: |
|
743 |
|
minimize c @ x |
|
744 |
|
such that: |
|
745 |
|
A_ub @ x <= b_ub |
|
746 |
|
A_eq @ x == b_eq |
|
747 |
|
lb <= x <= ub |
|
748 |
|
""" |
|
749 |
|
|
|
750 |
|
def __init__(self, points): |
|
751 |
|
self.gabriel = Gabriel(points) |
|
752 |
|
|
|
753 |
|
|
|
754 |
|
self.points = points |
|
755 |
|
nsim, ndim = points.shape |
|
756 |
|
size = nsim * (nsim - 1) // 2 |
|
757 |
|
self.condensed_contacts = np.zeros(size, dtype=bool) |
|
758 |
|
self.arange = np.arange(nsim) |
|
759 |
|
self.nsim = nsim |
|
760 |
|
|
|
761 |
|
self.A = A = np.empty((ndim + 1, nsim)) |
|
762 |
|
A[:-1] = points.T |
|
763 |
|
A[-1] = 1 |
|
764 |
|
|
|
765 |
|
self._b = np.empty(ndim + 1) |
|
766 |
|
self._b[-1] = 1 |
|
767 |
|
|
|
768 |
|
self.c = np.sum(np.square(points), axis=1) |
|
769 |
|
|
|
770 |
|
def LP(self, couple_indices): |
|
771 |
|
i, j = couple_indices |
|
772 |
|
X = self.points |
|
773 |
|
|
|
774 |
|
midpoint = np.mean(X[[i, j]], axis=0) |
|
775 |
|
b_eq = self._b |
|
776 |
|
b_eq[:-1] = midpoint |
|
777 |
|
|
|
778 |
|
return linprog(self.c, A_eq=self.A, b_eq=b_eq) |
|
779 |
|
|
|
780 |
|
def is_couple(self, couple_indices, tol=1e-5): |
|
781 |
|
i, j = couple_indices |
|
782 |
|
if self._get_contact(i, j): |
|
783 |
|
return True |
|
784 |
|
|
|
785 |
|
if self.gabriel.is_couple(couple_indices): |
|
786 |
|
self._set_contact(i, j) |
|
787 |
|
return True |
|
788 |
|
|
|
789 |
|
result = self.LP(couple_indices) |
|
790 |
|
mask = result.x > tol |
|
791 |
|
self._handle_polytope(self.arange[mask]) |
|
792 |
|
|
|
793 |
|
return mask[i] and mask[j] |
|
794 |
|
|
|
795 |
|
#X = self.points |
|
796 |
|
#test = (np.inner(X[i], X[i]) + np.inner(X[j], X[j])) / 2 |
|
797 |
|
#return test <= (result.fun + tol) |
|
798 |
|
|
|
799 |
|
|
|
800 |
|
|
|
801 |
|
def _handle_polytope(self, polytope): |
|
802 |
|
if len(polytope) > 2: |
|
803 |
|
i = polytope[0] |
|
804 |
|
jj = polytope[1:] |
|
805 |
|
self._handle_polytope(jj) |
|
806 |
|
for j in jj: |
|
807 |
|
self._set_contact(i, j) |
|
808 |
|
else: |
|
809 |
|
i, j = polytope |
|
810 |
|
self._set_contact(i, j) |
|
811 |
|
|
|
812 |
|
def _set_contact(self, x, y): |
|
813 |
|
i, j = min(x, y), max(x, y) |
|
814 |
|
entry = self.nsim * i + j - ((i + 2) * (i + 1)) // 2 |
|
815 |
|
self.condensed_contacts[entry] = True |
|
816 |
|
|
|
817 |
|
def _get_contact(self, x, y): |
|
818 |
|
i, j = min(x, y), max(x, y) |
|
819 |
|
entry = self.nsim * i + j - ((i + 2) * (i + 1)) // 2 |
|
820 |
|
return self.condensed_contacts[entry] |
815 |
821 |
|
|
|
822 |
|
|
|
823 |
|
|
|
824 |
|
|
|
825 |
|
class ConvexLinear(LinearSolver): |
|
826 |
|
def __init__(self, points, convex_solver=convex_sprite): |
|
827 |
|
super().__init__(points) |
816 |
828 |
|
|
|
829 |
|
nsim, ndim = points.shape |
|
830 |
|
self.lifted_points = np.empty((nsim, ndim + 1)) |
|
831 |
|
self.lifted_points[:, :ndim] = points |
|
832 |
|
# kind of datascience. feature space, quadratic kernel... |
|
833 |
|
self.lifted_points[:, -1] = np.sum(np.square(points), axis=1) |
|
834 |
|
|
|
835 |
|
self.convex_solver = convex_solver |
|
836 |
|
|
|
837 |
|
def is_couple(self, couple_indices, tol=1e-5, **kwargs): |
|
838 |
|
i, j = couple_indices |
|
839 |
|
if self._get_contact(i, j): |
|
840 |
|
return True |
|
841 |
|
|
|
842 |
|
if self.gabriel.is_couple(couple_indices): |
|
843 |
|
self._set_contact(i, j) |
|
844 |
|
return True |
|
845 |
|
|
|
846 |
|
if self.convex_solver(self.lifted_points, couple_indices, tol=tol, **kwargs): |
|
847 |
|
self._set_contact(i, j) |
|
848 |
|
return True |
|
849 |
|
|
|
850 |
|
result = self.LP(couple_indices) |
|
851 |
|
mask = result.x > tol |
|
852 |
|
self._handle_polytope(self.arange[mask]) |
|
853 |
|
|
|
854 |
|
return mask[i] and mask[j] |
|
855 |
|
|
|
856 |
|
|
|
857 |
|
|