-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathgeometry.py
More file actions
191 lines (161 loc) · 18.6 KB
/
geometry.py
File metadata and controls
191 lines (161 loc) · 18.6 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
from scipy.spatial import HalfspaceIntersection, ConvexHull, QhullError
import numpy as np
from itertools import accumulate, repeat, chain
from scipy.optimize import linprog
from debug_funcs import *
COMPARE_PRECISION = 2**(-20)
class D0RegionError(ValueError):
pass
class ExpectedQhullError(Exception):
pass
#Input: base_sig, sigs
#Returns: TL between base_sig and sigs as
# equations of the form (x_b - x)a + (z_b-z)c + (b_0(y_b-y) + (w_b - w)) = 0 as a matrix
# where (x_b, y_b, z_b, w_b) is the base_sig.
def compute_TLs(base_sig, sigs, b_val = 0):
diff = base_sig - sigs
#Note: tested speed comparison with vstack and array is ~25% faster
A = np.array([diff[:,0],diff[:,2],(b_val * diff[:,1] + diff[:,3])]).T
return A
def intersect_half_planes(base_point, halfspaces):
"""
Qhull seems to havve incorrect behavior when the base point is on the edge of the halfspaces.
- Speed can be increased if we can identify when this incorect behavior happens.
ex. hps = np.array([[1, -1, -830], [1, 0, -520], [0, 1, 310], [-1, 4, 1680], [1, -2, -1150], [0, -1, -320]])
base_point = np.array([512.85714, -317.14286])
(from seq: tmRNA_Neis.gono._TRW-242231_1-363)
- Fixed by the is close check in the try block.
"""
try:
if np.any(np.isclose(halfspaces[:,:-1].dot(base_point), -halfspaces[:,-1])):
raise(ExpectedQhullError)
HSI = HalfspaceIntersection(halfspaces, base_point, incremental=False)
except (ExpectedQhullError, QhullError):
norm_vector = np.reshape(np.linalg.norm(halfspaces[:, :-1], axis=1),(halfspaces.shape[0], 1))
c = np.zeros((halfspaces.shape[1],))
c[-1] = -1
A = np.hstack((halfspaces[:, :-1], norm_vector))
b = - halfspaces[:, -1:]
res = linprog(c, A_ub=A, b_ub=b, bounds=(None, None))
base_point = res.x[:-1]
radius = res.x[-1]
#Check if region is 0 area
if radius < COMPARE_PRECISION:
try:
endpoints, base_point = eval_1D_region(base_point, halfspaces)
return endpoints, base_point, True
except D0RegionError:
return None, base_point, True
else:
HSI = HalfspaceIntersection(halfspaces, base_point, incremental=False)
return HSI, base_point, False
"""
This determines the endpoints of a 1D region by first finding the parallel halfplanes that define it.
Let one of the halfplanes be (v, b_v) (v*p + b_v <= 0 for any point p in the halfplane)
We then solve the following LP:
max r
s.t.
s*p + r(v*s) <= b_s (for all halfplanes (s, b_s) which are not both parallel and adjacent to (v, b_v))
v*p = b_v (for the halfplane (v, b_v))
Then, given an optimal solution with point x, and radius r. We return,
the endpoints of our region:
x + v'r and x - v'r
"""
def eval_1D_region(base_point, halfspaces):
adjacent_idx = (np.abs(halfspaces[:,:-1].dot(base_point) + halfspaces[:,-1]) < COMPARE_PRECISION)
adjacent_hps = halfspaces[adjacent_idx]
#Find the indices of adjacent halfplanes that are parallel
slopes = {} #slope: index where slope is first found
parallel_idxs = [] #final indices of parallel halfplanes
parallel_val = None #slope of parallel halfplanes
for i, h in enumerate(adjacent_hps):
da, dc, _ = h
sign = int(np.prod(np.sign([da, dc])))
gcd = np.gcd(int(da), int(dc))
slope = (abs(int(da)//gcd), abs(int(dc)//gcd), sign)
#If we have already found 2 halfplanes which are parrallel, add any others (so that they aren't considered in the LP)
if (parallel_val != None) and (slope == parallel_val):
parallel_idxs.append(i)
continue
try:
#If we have already found an index with a the same slope as this index, save them and the slope value.
parallel_idxs = [slopes[slope], i]
parallel_val = slope
except KeyError:
#Otherwise add slope value.
slopes[slope] = i
#If there are no parallel halplanes, then region is 0 areas
if (parallel_val == None):
raise D0RegionError("No parallel halfspaces (region is not 1D) or base point is not in halfspaces.")
parallel_hps = adjacent_hps[parallel_idxs]
matching_hps = (halfspaces[:, None] == parallel_hps).all(-1).any(-1)
non_parallel_hps = halfspaces[~matching_hps]
c = np.zeros((non_parallel_hps.shape[1],)) #lp c
c[-1] = -1 #-1 linprog minimizes
v = parallel_hps[0][(0,1),]
A_eq = np.hstack((v[None,:], np.zeros((1,1)))) #lp equality eqs
b_eq = -parallel_hps[0][-1:] #lp equality vals
vt = parallel_hps[0][(1,0),] #vector for one of the parallel half planes (v, b_v)
vt[0] = -vt[0]
dot_vector = np.abs(non_parallel_hps[:, :-1].dot(vt))
A_ub = np.hstack((non_parallel_hps[:, :-1], dot_vector[:,None])) #lp upper bound eqs
b_ub = -non_parallel_hps[:,-1:] #lp upper bound vals
res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=(None, None))
base_point = res.x[:-1]
radius = res.x[-1]
if radius < COMPARE_PRECISION:
raise D0RegionError(f"Region is a point, not 1D, {base_point}")
endpoints = (base_point + radius * vt), (base_point - radius * vt)
# nice_print_points(endpoints, hull=False)
return endpoints, base_point
if __name__ == "__main__":
# bs = np.array((0, 0, 0, -2980))
# bp = np.array((1036,1871))
# # np.random.seed(1)
# # sigs = np.random.randint(0,3,size = (10,4))
# sigs = np.array(list(set([(1, 7, 3, -4010),(1, 8, 4, -3980),(1, 8, 4, -3980),(1, 4, 5, -3490),(1, 3, 7, -2430),(1, 4, 5, -3490),(1, 7, 3, -4010),(2, 6, 6, -4070),(1, 8, 4, -3980),(2, 6, 6, -4070),(3, 3, 9, -3970),(1, 8, 4, -3980),(3, 3, 9, -3970),(4, 3, 12, -3430),(3, 4, 10, -3370),(3, 3, 9, -3970),(3, 4, 10, -3370),(2, 4, 7, -3880),(3, 3, 9, -3970),(1, 8, 4, -3980),(2, 4, 7, -3880),(4, 3, 12, -3430),(5, 2, 15, -2660),(4, 2, 13, -2690),(4, 3, 12, -3430),(3, 4, 10, -3370),(4, 2, 13, -2690),(5, 2, 15, -2660),(6, 1, 18, -1690),(6, 0, 19, -800),(5, 2, 15, -2660),(4, 2, 13, -2690),(4, 2, 14, -1840),(6, 1, 18, -1690),(7, 0, 21, -480),(6, 0, 19, -800),(7, 0, 21, -480),(8, 4, 24, 940),(7, 1, 22, 440),(7, 0, 21, -480),(6, 0, 19, -800),(7, 1, 22, 440),(8, 4, 24, 940),(9, 2, 27, 2650),(8, 3, 25, 2050),(8, 4, 24, 940),(7, 1, 22, 440),(8, 3, 25, 2050),(9, 2, 27, 2650),(7, 1, 24, 2940),(8, 3, 25, 2050),(7, 1, 24, 2940),(8, 3, 25, 2050),(6, 0, 22, 2420),(7, 1, 24, 2940),(5, 0, 21, 3630),(6, 0, 22, 2420),(5, 0, 21, 3630),(3, 0, 17, 3430),(6, 0, 22, 2420),(3, 0, 17, 3430),(6, 0, 22, 2420),(3, 0, 16, 2130),(3, 0, 17, 3430),(2, 1, 15, 3530),(1, 2, 12, 2260),(3, 0, 17, 3430),(1, 2, 12, 2260),(2, 1, 14, 2180),(3, 0, 17, 3430),(3, 0, 16, 2130),(2, 1, 14, 2180),(2, 1, 15, 3530),(1, 2, 13, 3720),(1, 2, 12, 2260),(1, 2, 12, 2260),(2, 1, 14, 2180),(2, 3, 13, 1120),(1, 2, 12, 2260),(1, 4, 10, 190),(2, 2, 11, -990),(1, 2, 12, 2260),(2, 2, 11, -990),(2, 2, 12, 60),(1, 4, 10, 190),(1, 3, 7, -2430),(2, 2, 11, -990),(1, 3, 7, -2430),(3, 4, 10, -3370),(2, 4, 7, -3880),(1, 3, 7, -2430),(3, 4, 10, -3370),(2, 3, 10, -1800),(1, 3, 7, -2430),(2, 2, 11, -990),(2, 3, 10, -1800),(1, 3, 7, -2430),(1, 4, 5, -3490),(2, 4, 7, -3880),(1, 8, 4, -3980),(1, 4, 5, -3490),(2, 4, 7, -3880),(3, 4, 10, -3370),(4, 2, 13, -2690),(2, 3, 10, -1800),(4, 2, 13, -2690),(4, 2, 14, -1840),(3, 2, 12, -1860),(4, 2, 13, -2690),(2, 3, 10, -1800),(3, 2, 12, -1860),(6, 0, 19, -800),(7, 1, 22, 440),(6, 0, 22, 2420),(6, 0, 19, -800),(4, 1, 16, 70),(5, 1, 18, 160),(6, 0, 19, -800),(3, 1, 13, -990),(4, 1, 16, 70),(7, 1, 22, 440),(8, 3, 25, 2050),(6, 0, 22, 2420),(6, 0, 22, 2420),(4, 1, 17, 1140),(5, 1, 18, 160),(2, 2, 11, -990),(2, 3, 10, -1800),(3, 2, 12, -1860),(2, 2, 11, -990),(3, 2, 12, -1860),(3, 1, 13, -990),(2, 2, 11, -990),(3, 1, 13, -990),(3, 1, 14, 40),(2, 2, 11, -990),(3, 1, 14, 40),(3, 2, 15, 1070),(2, 2, 11, -990),(2, 2, 12, 60),(3, 2, 15, 1070),(4, 2, 14, -1840),(3, 2, 12, -1860),(3, 1, 13, -990),(3, 0, 16, 2130),(4, 1, 17, 1140),(4, 1, 16, 70),(3, 0, 16, 2130),(4, 1, 16, 70),(3, 2, 15, 1070),(3, 1, 13, -990),(4, 1, 16, 70),(3, 1, 14, 40),(4, 1, 17, 1140),(4, 1, 16, 70),(5, 1, 18, 160),(2, 3, 13, 1120),(2, 2, 12, 60),(3, 2, 15, 1070),(4, 1, 16, 70),(3, 1, 14, 40),(3, 2, 15, 1070)])))
# transform = np.array([[1,0,0,0],[0,1,0,0],[-3,0,1,0],[0,0,0,1]])
# sigs = sigs @ transform.T
# # print(sigs)
# halfspaces = compute_TLs(bs, sigs)
# halfspaces = np.vstack([halfspaces, np.array([1,1,-3000])])
# region, HSI = intersect_half_planes(bp, halfspaces)
# print(HSI.intersections)
# contracted = contract_hull(region, HSI.intersections)
"""
1.0x + -1.0y + -830.0 <= 0
1.0x + 0.0y + -520.0 <= 0
0.0x + 1.0y + 310.0 <= 0
-1.0x + 4.0y + 1680.0 <= 0
1.0x + -2.0y + -1150.0 <= 0
0.0x + -1.0y + -320.0 <= 0
"""
hps = np.array([[26.0, 28.0, 21620.0],[16.0, -30.0, -17010.0],[25.0, -39.0, -18300.0],[8.0, 18.0, 21930.0],[-30.0, 28.0, -9010.0],[-7.0, 13.0, 5720.0],[17.0, -10.0, 6080.0],[25.0, -6.0, 11200.0],[20.0, 19.0, 24810.0],[20.0, -34.0, -17310.0],[18.0, -32.0, -17090.0],[19.0, -33.0, -17190.0],[17.0, 6.0, 18000.0],[21.0, 9.0, 20220.0],[12.0, 22.0, 25560.0],[-7.0, 28.0, 18460.0],[21.0, -3.0, 13230.0],[14.0, 28.0, 27650.0],[17.0, -31.0, -17020.0],[4.0, 11.0, 14310.0],[9.0, -15.0, -6290.0],[12.0, 7.0, 16750.0],[-15.0, 9.0, -8380.0],[21.0, 26.0, 25490.0],[20.0, 28.0, 25680.0],[17.0, 27.0, 27000.0],[16.0, 28.0, 27180.0],[15.0, 28.0, 27430.0],[18.0, -28.0, -11320.0],[1.0, -12.0, -12490.0],[10.0, -20.0, -10950.0],[14.0, -23.0, -9840.0],[11.0, -20.0, -9710.0],[-10.0, -3.0, -18320.0],[-24.0, 17.0, -13860.0],[17.0, -30.0, -15040.0],[-1.0, -12.0, -15810.0],[-13.0, 2.0, -14840.0],[-7.0, -4.0, -13400.0],[-18.0, 8.0, -15350.0],[-28.0, 24.0, -10080.0],[-21.0, 15.0, -9950.0],[-13.0, 1.0, -16990.0],[23.0, -37.0, -17820.0],[20.0, -33.0, -15440.0],[23.0, -36.0, -16060.0],[19.0, -30.0, -12680.0],[17.0, 17.0, 24100.0],[16.0, 24.0, 27040.0],[17.0, 24.0, 26820.0],[19.0, 27.0, 26290.0],[18.0, 28.0, 26500.0],[-22.0, 15.0, -12030.0],[-27.0, 22.0, -11380.0],[23.0, 27.0, 24130.0],[22.0, -26.0, -5630.0],[17.0, -20.0, -3630.0],[18.0, 22.0, 26080.0],[17.0, 23.0, 26570.0],[20.0, 21.0, 25310.0],[20.0, 17.0, 24110.0],[18.0, 23.0, 26330.0],[21.0, 23.0, 25210.0],[22.0, 18.0, 23710.0],[25.0, 23.0, 22450.0],[25.0, 9.0, 18510.0],[24.0, 2.0, 16100.0],[23.0, 0.0, 15030.0],[23.0, -5.0, 12030.0],[19.0, -27.0, -9180.0],[24.0, -20.0, 1240.0],[22.0, -20.0, 260.0],[18.0, -18.0, -720.0],[20.0, -18.0, 790.0],[18.0, -14.0, 3100.0],[18.0, -12.0, 4960.0],[20.0, -16.0, 2630.0],[19.0, -24.0, -5970.0],[1.0, 20.0, 18970.0],[16.0, -12.0, 3400.0],[21.0, -10.0, 8350.0],[18.0, -1.0, 13830.0],[22.0, 7.0, 19050.0],[19.0, 23.0, 25980.0],[20.0, 22.0, 25460.0],[25.0, 16.0, 20820.0],[25.0, 11.0, 19230.0],[23.0, 12.0, 21050.0],[24.0, 12.0, 20480.0],[23.0, 9.0, 19750.0],[21.0, 16.0, 23560.0],[25.0, 10.0, 18880.0],[25.0, 12.0, 19550.0],[24.0, 15.0, 21450.0],[22.0, -15.0, 4660.0],[20.0, -12.0, 6150.0],[17.0, -6.0, 9480.0],[9.0, 8.0, 15750.0],[6.0, 16.0, 19600.0],[9.0, 14.0, 19900.0],[12.0, 15.0, 21810.0],[14.0, 15.0, 22470.0],[11.0, 9.0, 17590.0],[11.0, 14.0, 20800.0],[3.0, 16.0, 17510.0],[-3.0, 21.0, 16580.0],[-18.0, 28.0, 8110.0],[-26.0, 27.0, -2820.0],[-16.0, 16.0, -1460.0],[-3.0, 0.0, -3380.0],[0.0, -8.0, -8850.0],[7.0, -16.0, -9730.0],[-9.0, 5.0, -5060.0],[-3.0, -3.0, -6750.0],[-30.0, 27.0, -10660.0],[-7.0, -1.0, -9530.0],[14.0, 2.0, 14130.0],[17.0, -3.0, 11870.0],[-6.0, -8.0, -19600.0],[-3.0, -9.0, -14340.0],[4.0, 26.0, 24550.0],[5.0, 28.0, 25870.0],[-2.0, 28.0, 22070.0],[-5.0, 28.0, 20000.0],[-4.0, 25.0, 18810.0],[-2.0, 26.0, 20930.0],[-3.0, 28.0, 21440.0],[-4.0, 27.0, 20160.0],[2.0, 28.0, 24340.0],[0.0, 28.0, 23210.0],[-7.0, 26.0, 17050.0],[-12.0, 23.0, 9880.0],[13.0, -24.0, -12230.0],[15.0, -26.0, -12320.0],[14.0, -25.0, -12270.0],[12.0, -25.0, -14850.0],[24.0, -34.0, -12780.0],[21.0, -11.0, 7550.0],[18.0, -7.0, 9230.0],[12.0, 18.0, 23560.0],[10.0, 17.0, 22210.0],[4.0, -2.0, 2170.0],[23.0, 28.0, 23870.0],[20.0, 25.0, 25810.0],[19.0, -21.0, -2910.0],[22.0, -23.0, -2600.0],[-6.0, 28.0, 19260.0],[12.0, -23.0, -12210.0],[20.0, -21.0, -2130.0],[21.0, -32.0, -13020.0],[21.0, -5.0, 11940.0],[21.0, -9.0, 9100.0],[23.0, -13.0, 6630.0],[23.0, -9.0, 9470.0],[25.0, -16.0, 4530.0],[25.0, -12.0, 7340.0],[4.0, -8.0, -4120.0],[7.0, -13.0, -6270.0],[16.0, 15.0, 22920.0],[12.0, 1.0, 12070.0],[-29.0, 26.0, -9460.0],[23.0, 26.0, 24320.0],[23.0, 24.0, 24180.0],[24.0, 24.0, 23410.0],[19.0, -12.0, 5580.0],[21.0, -34.0, -15610.0],[0.0, -4.0, -4300.0],[-3.0, -6.0, -10270.0],[-5.0, -2.0, -8110.0],[-3.0, -4.0, -7920.0],[21.0, 11.0, 21190.0],[19.0, 16.0, 23690.0],[20.0, 12.0, 21720.0],[-18.0, 10.0, -11910.0],[25.0, -29.0, -6840.0],[-5.0, -5.0, -11800.0],[-2.0, -4.0, -6670.0],[24.0, -13.0, 6720.0],[25.0, -15.0, 5260.0],[25.0, -13.0, 6660.0],[21.0, 18.0, 24280.0],[21.0, 17.0, 23920.0],[20.0, 26.0, 25890.0],[-5.0, 9.0, 3900.0],[-1.0, 10.0, 8900.0],[21.0, -29.0, -9630.0],[21.0, -31.0, -11880.0],[24.0, -36.0, -15050.0],[25.0, -37.0, -15390.0],[24.0, -37.0, -16360.0],[-28.0, 28.0, -4860.0],[-24.0, 28.0, 940.0],[9.0, 23.0, 25210.0],[9.0, 28.0, 27110.0],[-6.0, 3.0, -3610.0],[-9.0, 10.0, 420.0],[6.0, 15.0, 18920.0],[4.0, 16.0, 18250.0],[5.0, 16.0, 18960.0],[23.0, -10.0, 8790.0],[22.0, -7.0, 10710.0],[16.0, 1.0, 14410.0],[17.0, 5.0, 17360.0],[18.0, 0.0, 14480.0],[22.0, 28.0, 24530.0],[22.0, 27.0, 24730.0],[21.0, 28.0, 25110.0],[-2.0, 8.0, 6020.0],[-3.0, 12.0, 8820.0],[-9.0, -3.0, -15400.0],[10.0, -8.0, 1970.0],[13.0, -17.0, -4340.0],[5.0, 2.0, 7060.0],[0.0, 5.0, 5130.0],[4.0, 3.0, 7080.0],[0.0, 8.0, 7970.0],[15.0, -22.0, -7660.0],[11.0, -17.0, -6350.0],[17.0, -17.0, -580.0],[-6.0, -6.0, -14740.0],[-4.0, -8.0, -14460.0],[-4.0, -7.0, -12980.0],[15.0, -8.0, 6320.0],[8.0, -1.0, 6910.0],[9.0, -6.0, 3020.0],[15.0, -9.0, 5390.0],[19.0, -7.0, 9750.0],[18.0, -6.0, 10040.0],[6.0, -2.0, 4110.0],[6.0, 0.0, 6070.0],[0.0, -9.0, -10020.0],[-2.0, -7.0, -10180.0],[-15.0, 5.0, -14000.0],[-17.0, 7.0, -14900.0],[-8.0, 4.0, -4940.0],[-6.0, -1.0, -8230.0],[-9.0, 6.0, -3900.0],[9.0, -17.0, -8550.0],[11.0, 12.0, 19550.0],[20.0, 7.0, 19040.0],[12.0, -20.0, -8620.0],[2.0, 25.0, 22880.0],[2.0, 26.0, 23410.0],[2.0, 23.0, 21640.0],[4.0, 23.0, 22860.0],[7.0, 23.0, 24380.0],[8.0, 22.0, 24260.0],[7.0, 17.0, 20830.0],[7.0, 19.0, 22050.0],[4.0, 22.0, 22230.0],[5.0, 17.0, 19640.0],[3.0, 22.0, 21630.0],[3.0, 19.0, 19660.0],[4.0, 17.0, 18970.0],[23.0, -8.0, 10130.0],[24.0, -12.0, 7420.0],[25.0, 0.0, 14670.0],[24.0, -3.0, 13210.0],[25.0, 5.0, 16900.0],[24.0, 9.0, 19280.0],[25.0, 6.0, 17310.0],[22.0, 19.0, 23960.0],[21.0, 19.0, 24530.0],[23.0, 17.0, 22710.0],[-3.0, -10.0, -16090.0],[23.0, 6.0, 18330.0],[24.0, -16.0, 4480.0],[24.0, -14.0, 6010.0],[25.0, -17.0, 3800.0],[25.0, -22.0, -260.0],[15.0, 17.0, 23720.0],[15.0, -14.0, 580.0],[12.0, -10.0, 1810.0],[16.0, -17.0, -1480.0],[-5.0, -8.0, -16450.0],[-9.0, -4.0, -17680.0],[8.0, 28.0, 26870.0],[-26.0, 28.0, -1870.0],[4.0, 10.0, 13460.0],[-27.0, 25.0, -6800.0],[-24.0, 23.0, -4470.0],[-28.0, 25.0, -8560.0],[-29.0, 27.0, -8040.0],[-28.0, 26.0, -7220.0],[-29.0, 28.0, -6750.0],[-21.0, 28.0, 4700.0],[-19.0, 28.0, 6980.0],[-19.0, 24.0, 3310.0],[-22.0, 24.0, -490.0],[-19.0, 26.0, 5190.0],[-20.0, 26.0, 4030.0],[-15.0, 22.0, 5830.0],[-13.0, 18.0, 4130.0],[-15.0, 17.0, 790.0],[-15.0, 23.0, 6760.0],[-15.0, 16.0, -260.0],[-11.0, 13.0, 1260.0],[-9.0, 12.0, 2520.0],[-10.0, 13.0, 2420.0],[-19.0, 16.0, -5530.0],[-20.0, 16.0, -7010.0],[0.0, 26.0, 22210.0],[-2.0, 27.0, 21520.0],[-3.0, 27.0, 20850.0],[11.0, -8.0, 2860.0],[20.0, 5.0, 17900.0],[21.0, 7.0, 19080.0],[22.0, 1.0, 15600.0],[13.0, 28.0, 27800.0],[14.0, 27.0, 27660.0],[12.0, -24.0, -13520.0],[14.0, -26.0, -13580.0],[-1.0, 1.0, -40.0],[2.0, -2.0, 80.0],[2.0, -4.0, -2030.0],[3.0, -7.0, -4160.0],[5.0, -11.0, -6330.0],[1.0, -3.0, -2070.0],[24.0, 23.0, 23330.0],[25.0, 24.0, 22530.0],[24.0, 19.0, 22430.0],[24.0, 16.0, 21700.0],[5.0, -18.0, -15260.0],[14.0, 8.0, 18230.0],[14.0, 9.0, 18870.0],[25.0, -33.0, -10920.0],[-18.0, 11.0, -10430.0],[5.0, 27.0, 25500.0],[5.0, -10.0, -5200.0],[4.0, -9.0, -5240.0],[7.0, 26.0, 25940.0],[6.0, 27.0, 26000.0],[5.0, 26.0, 25050.0],[10.0, 15.0, 20980.0],[25.0, -21.0, 620.0],[24.0, -21.0, 360.0],[-24.0, 27.0, -20.0],[-22.0, 27.0, 2540.0],[-22.0, 28.0, 3500.0],[-21.0, 27.0, 3770.0],[-23.0, 27.0, 1270.0],[24.0, -28.0, -6400.0],[22.0, -24.0, -3590.0],[13.0, 25.0, 26960.0],[14.0, 26.0, 27400.0],[21.0, -15.0, 4150.0],[20.0, -14.0, 4440.0],[21.0, -13.0, 5860.0],[22.0, -13.0, 6370.0],[23.0, -16.0, 4290.0],[13.0, 16.0, 22760.0],[0.0, -11.0, -12570.0],[-5.0, 25.0, 18000.0],[22.0, -6.0, 11370.0],[20.0, 20.0, 25060.0],[18.0, -2.0, 13110.0],[17.0, -1.0, 13400.0],[19.0, -2.0, 13470.0],[19.0, 0.0, 14760.0],[18.0, 1.0, 15120.0],[23.0, 16.0, 22450.0],[6.0, -17.0, -12270.0],[3.0, -14.0, -12390.0],[11.0, -24.0, -14870.0],[7.0, -19.0, -13620.0],[9.0, -22.0, -14970.0],[-5.0, 27.0, 19420.0],[-6.0, 26.0, 17920.0],[-4.0, -9.0, -16260.0],[22.0, -11.0, 7930.0],[22.0, -10.0, 8640.0],[-2.0, 10.0, 7920.0],[8.0, -22.0, -17250.0],[-2.0, -12.0, -18430.0],[-4.0, -10.0, -18890.0],[12.0, -26.0, -17090.0],[10.0, -24.0, -17150.0],[17.0, 28.0, 26850.0],[14.0, 17.0, 23500.0],[24.0, 0.0, 14970.0],[25.0, -3.0, 13000.0],[25.0, -4.0, 12410.0],[20.0, -30.0, -11640.0],[20.0, 24.0, 25710.0],[9.0, -2.0, 6830.0],[23.0, -30.0, -9140.0],[24.0, -33.0, -11660.0],[18.0, 27.0, 26650.0],[-23.0, 20.0, -6640.0],[24.0, -15.0, 5260.0],[25.0, -18.0, 3020.0],[25.0, -38.0, -16700.0],[13.0, 9.0, 18470.0],[12.0, 10.0, 18690.0],[-4.0, 1.0, -3450.0],[-5.0, 7.0, 1830.0],[-6.0, 8.0, 1760.0],[-7.0, 9.0, 1680.0],[-8.0, 10.0, 1580.0],[0.0, 1.0, -10000.0],[-1.0, 0.0, -10000.0],[0.0, -1.0, -10000.0],[1.0, 0.0, -10000.0]])
base_point = np.array([-1066, -1026])
hps, bp, zero = intersect_half_planes(base_point, hps)
print(hps, bp, zero)
# fig = plt.figure()
# ax = fig.add_subplot(1, 1, 1) #aspect='equal'
# fmt = {"color": None, "edgecolor": "b", "alpha": 0.5}
# x = np.linspace(-3000,3000)
# for h in halfspaces:
# hlist = h.tolist()
# if h[1]== 0:
# ax.axvline(-h[2]/h[0], label='{}x+{}y+{}=0'.format(*hlist),color="black")
# # xi = np.linspace(xlim[sign], -h[2]/h[0], 100)
# # ax.fill_between(xi, ylim[0], ylim[1], **fmt)
# else:
# ax.plot(x, (-h[2]-h[0]*x)/h[1], label='{}x+{}y+{}=0'.format(*hlist),color="black")
# # ax.fill_between(x, (-h[2]-h[0]*x)/h[1], ylim[sign], **fmt)
# # x, y = zip(*halfspaces.intersections)
# for h in region:
# hlist = h.tolist()
# if h[1]== 0:
# ax.axvline(-h[2]/h[0], label='{}x+{}y+{}=0'.format(*hlist),color="red")
# # xi = np.linspace(xlim[sign], -h[2]/h[0], 100)
# # ax.fill_between(xi, ylim[0], ylim[1], **fmt)
# else:
# ax.plot(x, (-h[2]-h[0]*x)/h[1], label='{}x+{}y+{}=0'.format(*hlist),color="red")
# # ax.fill_between(x, (-h[2]-h[0]*x)/h[1], ylim[sign], **fmt)
# # x, y = zip(*halfspaces.intersections)
# plt.show()