PYTHON   92

# LP

Guest on 14th June 2022 01:30:08 AM

1.
2. from frac import *
3. from frac_matrix import *
4. from array_strings import *
5.
6. import sys
7.
8. Infinity = None
9.
10. class LP_Dictionary:
11.
12.         # subscripts index rows of F, indices its columns
13.         """ Data:
14.
15.                 An M x N+1 matrix F of affine functions with constant in last place
16.                         These correspond to inequalities
17.                 An N+1 x N+1 affine transformation matrix R: original F.R = current F
18.                 A list of coordinate functions
19.                 A list of other functions in the inequalities
20.                 An index, listing for each basis variable which column it corresponds to
21.         """
22.         # args = R, basis, others, index
23.         def __init__(self, F, *args):
24.                 M = len(F)
25.                 if M == 0:
26.                         print "LP_Dictionary: Matrix must not be empty!"
27.                         sys.exit(1)
28.                 if len(args) > 0:
29.                         # we already have a vertex frame
30.                         self.F = F
31.                         self.R = args
32.                         self.basis = args
33.                         self.others = args
34.                         self.index = args
35.                         self.dim = None
36.                         self.rank = len(F) - 1
37.                 else:
38.                         N = len(F)-1 # number of variables
39.                         # tack R onto the bottom of F
40.                         R = matrix.identity(N+1)
41.                         for x in R:
42.                                 F.append(x)
43.                         [F, self.index, self.basis, self.others] = self.gauss(F, M, N)
44.                         self.original_R = F[M:]
45.                         F = F[:M]
46.                         rank = len(self.basis)
47.                         # delete penultimate N-rank rows
48.                         for i in range(M):
49.                                 f = F[i]
50.                                 F[i] = f[:rank] + f[-1:]
51.                         self.F = F
52.                         self.R = matrix.identity(rank+1)
53.                         self.rank = rank
54.                         self.null_dim = N - rank
55.                         # immediately move to vertex frame if possible
56.                         X = self.to_vertex_frame()
57.                         if X is None:
58.                                 self.dim = -1
59.                         else:
60.                                 [ self.F, self.R, self.basis, self.others, self.index ] = X
61.                                 self.dim = None
62.
63.         def clone(self):
64.                 F = []
65.                 for x in self.F:
66.                         f = []
67.                         for y in x:
68.                                 f.append(y)
69.                         F.append(f)
70.                 R = []
71.                 for x in self.R:
72.                         f = []
73.                         for y in x:
74.                                 f.append(y)
75.                         R.append(f)
76.                 I = []
77.                 for x in self.index:
78.                         I.append(x)
79.                 B = []
80.                 for x in self.basis:
81.                         B.append(x)
82.                 O = []
83.                 for x in self.others:
84.                         O.append(x)
85.                 return LP_Dictionary(F, R, B, O, I)
86.
87.         def pivot(self, e, E):
88.                 # i.e we are swapping E and e
89.                 # e is one of the basis variables
90.                 # this follows along the e-axis to the next vertex
91.                 # and we have chosen a variable E to replace it in the basis
93.                 [F, R, basis, others, index] = [ self.F, self.R, self.basis, self.others, self.index]
94.                 # E replaces e as basis variable
95.                 k = index[e]
96.                 index[e] = None
97.                 basis[k] = E
98.                 index[E] = k
99.                 # the column of E in F is k
100.                 others.append(e)
101.                 others.remove(E) # we don't care about order in others
102.
103.                 M = len(F)
104.                 N = len(F)-1
105.                 p = F[E][k]
106.                 for i in range(M):
107.                         F[i][k] /= p
108.                 for i in range(N):
109.                         R[i][k] /= p
110.                 # subtract a column[a][m] * column[e] from column[a], all a != e
111.                 for a in range(k):
112.                         p = F[E][a]
113.                         for b in range(M):
114.                                 F[b][a] -= p*F[b][k]
115.                         for b in range(N):
116.                                 R[b][a] -= p*R[b][k]
117.                 for a in range(k+1, N+1):
118.                         p = F[E][a]
119.                         for b in range(M):
120.                                 F[b][a] -= p*F[b][k]
121.                         for b in range(N):
122.                                 R[b][a] -= p*R[b][k]
123.
124.         # have to use technique in the notes because of possibility of empty set
125.         def to_vertex_frame(self):
126.                 F = self.F
127.                 M = len(F)
128.                 N = len(F) - 1
129.                 basis = self.basis
130.                 # make a new copy of F and a new polyhedron
131.                 # add two functions x_0, -x_0 + 1 in notation of notes
132.                 # homogenize current basis
133.                 H = []
134.                 for x in F:
135.                         row = []
136.                         for y in x:
137.                                 row.append(y)
138.                         row.append(frac(0))
139.                         H.append(row)
140.                 H.append([frac(0)]*(N+2))
141.                 H[-1][N] = frac(1) # x_0
142.                 H.append([frac(0)]*(N+2))
143.                 H[-1][N] = frac(-1)
144.                 H[-1][N+1] = frac(1) # -x_0 + 1
145.                 I = self.index + [len(basis), None]
146.                 B = self.basis + [M]
147.                 O = self.others + [M+1]
148.                 R = matrix.identity(N+2)
149.                 X = LP_Dictionary(H, R, B, O, I)
150.                 z = [ x for x in H[M] ]
151.                 X.to_maximal_vertex(z)
152.                 # either 0 or 1
153.                 # if 1: we are succesful
154.                 # if 0: polyhedron is empty
155.                 c = H[M][-1]
156.                 if c == 0:
157.                         return None # this system is empty
158.
159.                 m = X.index[M+1]
160.                 # now chop H to get the vertex frame
161.                 H = H[:-3]
162.                 I = []
163.                 B = []
164.                 for x in X.basis:
165.                         if x != M+1:
166.                                 B.append(x)
167.                 for x in X.index[:-1]:
168.                         if x < m:
169.                                 I.append(x)
170.                         else:
171.                                 I.append(x-1)
172.                 nF = []
173.                 for i in range(len(H)):
174.                         row = []
175.                         for j in range(len(H)):
176.                                 if j != m:
177.                                         row.append(H[i][j])
178.                         nF.append(row)
179.                 nR = []
180.                 for i in range(len(X.R)):
181.                         row =[]
182.                         if i != len(X.R)-2:
183.                                 nR.append(row)
184.                                 for j in range(len(X.R)):
185.                                         if j != m:
186.                                                 row.append(X.R[i][j])
187.
188.                 O = range(M)
189.                 for x in B:
190.                         O.remove(x)
191.                 self.others = O
192.                 return [ nF, matrix(nR), B, O, I ]
193.
194.
195.         # finds first that increases f
196.         def exiting_variable(self, f):
197.                 basis = sorted(self.basis)
198.                 for x in basis:
199.                         i = self.index[x]
200.                         if f[i] > 0:
201.                                 return x
202.                 return None # we are at a maximum
203.
204.         def entering_variable(self, x):
205.                 others = sorted(self.others)
206.                 F = self.F
207.                 N = len(F) - 1
208.                 m = None
209.                 # find minimum value
210.                 minimum = None
211.                 for y in others:
212.                         phi = F[y]
213.                         i = self.index[x]
214.                         if phi[i] < 0:
215.                                 c = -phi[N]/phi[i]
216.                                 if minimum is not None and c < minimum:
217.                                         m = y
218.                                         minimum = c
219.                                 elif minimum is None:
220.                                         m = y
221.                                         minimum = c
222.                 if minimum is None:
223.                         # system unbounded, the edge x goes to infinity
224.                         return None
225.                 return m
226.
227.         def to_maximal_vertex(self, z):
228.                 # z in which coordinates?
229.                 self.F.append(z)
230.                 x = self.exiting_variable(z)
231.                 while x is not None:
232.                         E = self.entering_variable(x) # so E is smallest with F[E][e] < 0 minimal
233.                         if E is None:
234.                                 return Infinity
235.                         self.pivot(x, E)
236.                         x = self.exiting_variable(z)
237.                 return z[-1]
238.
239.         def maximum_value(self, z):
240.                 X = self.clone()
241.                 # clone self
242.                 X.F.append(z)
243.                 x = X.exiting_variable(z)
244.                 while x is not None:
245.                         E = X.entering_variable(x) # so E is smallest with F[E][e] < 0 minimal
246.                         if E is None:
247.                                 return Infinity
248.                         X.pivot(x, E)
249.                         x = X.exiting_variable(z)
250.                 return z[-1]
251.
252.         def tangent_cone(self):
253.                 H = []
254.                 I = []
255.                 B = []
256.                 n = 0
257.                 for x in self.basis:
258.                         H.append(self.F[x])
259.                         I.append(self.index[x])
260.                         B.append(n)
261.                         n += 1
262.                 rank = n
263.                 O = []
264.                 for x in self.others:
265.                         f = self.F[x]
266.                         if f[-1] == 0:
267.                                 H.append(self.F[x])
268.                                 I.append(None)
269.                                 O.append(n)
270.                                 n += 1
271.                 R = matrix(rank+1, rank+1)
272.                 for i in range(rank):
273.                         R[i][i] = frac(1)
274.                 cone = LP_Dictionary(H, R, B, O, I)
275.                 return cone
276.
277.         def dimension(self):
278.                 if self.dim is not None:
279.                         return self.dim
280.                 return self._dimension_(self.tangent_cone())
281.
282.         # or static?, so easy to call recursively?
283.         def _dimension_(cone):
284.                 ## print "-------------------------------"
285.                 ## cone._print_()
286.                 # assumed to be a vertex frame
287.                 rank = len(cone.basis)
288.                 if rank == 1:
289.                         k = cone.basis[-1]
290.                         f = vector(cone.F[k]).clone()
291.                         C = cone.maximum_value(f)
292.                         if C is Infinity:
293.                                 return 1
294.                         else: #C = 0
295.                                 return 0
296.                 k = cone.basis[-1]
297.                 f = vector(cone.F[k]).clone()
298.                 C = cone.maximum_value(f)
299.                 if C == 0:
300.                         # cut down to one dimension one less
301.                         # make up a new cone d by setting x_{m} = 0 in all functions
302.                         # discarding x_{m} as a variable
303.                         m = cone.index[k] # f_{k} = x_{m}
304.                         H = []
305.                         # each row of F except k goes into H
306.                         # and we remove column m
307.                         for i in range(len(cone.F)):
308.                                 if i != k:
309.                                         row = []
310.                                         f = cone.F[i]
311.                                         for j in range(len(f)):
312.                                                 if j != m:
313.                                                         row.append(f[j])
314.                                         H.append(row)
315.                         N =len(H)-1
316.                         R = matrix(N, N+1)
317.                         B = cone.basis[:-1]
318.                         I = []
319.                         for i in range(len(cone.index)):
320.                                 if i != k:
321.                                         x = cone.index[i]
322.                                         if x is not None:
323.                                                 if x < m:
324.                                                         I.append(x)
325.                                                 else:
326.                                                         I.append(x-1)
327.                                         else:
328.                                                 I.append(x)
329.                         O = range(len(H))
330.                         for x in B:
331.                                 O.remove(x)
332.                         D = LP_Dictionary(H, R, B, O, I)
333.                         return D.dimension()
334.                 else: # C = Infinity
335.                         # make a new cone D of dimension one less
336.                         # by discarding all functions that depend on x_{m}
337.                         m = cone.index[k] # f_{k} = x_{m}
338.                         H = []
339.                         n = 0
340.                         tr = [None]*len(cone.F) # translation, old to new
341.                         for i in range(len(cone.F)):
342.                                 f = cone.F[i]
343.                                 c = f[m]
344.                                 if c == 0:
345.                                         tr[i] = n
346.                                         n += 1
347.                                         h = []
348.                                         for j in range(len(f)):
349.                                                 if j != m:
350.                                                         h.append(f[j])
351.                                         H.append(h)
352.                         M = len(H)
353.                         N = len(H) - 1
354.                         B = []
355.                         for x in cone.basis:
356.                                 if x != k:
357.                                         n = tr[x]
358.                                         B.append(n)
359.                         I = [None]*M
360.                         for i in range(len(B)):
361.                                 n = tr[i]
362.                                 x = cone.index[n]
363.                                 if x < m:
364.                                         I[n] = x
365.                                 else:
366.                                         I[n] = x-1
367.
368.                         R = matrix.identity(N+1)
369.                         O = range(M)
370.                         for x in B:
371.                                 O.remove(x)
372.                         D = LP_Dictionary(H, R, B, O, I)
373.
374.                         return 1 + D.dimension()
375.         _dimension_=staticmethod(_dimension_)
376.
377.         def _print_(self):
378.                 print_matrix(self.F)
379.                 print_matrix(self.R)
380.                 print self.basis, self.others, self.index
381.
382.         def __str__(self):
383.                 s = ""
384.                 return s
385.
386.
387.         # -----------------------------------------------------------------------
388.
389.         """ Performs gauss elimination on F,
390.         but by looking only in the upper left M x N for linearly independent rows
391.         but pivoting on the entire matrix.
392.         if lower (N+1) x (N+1) starts as I, it will end as affine R such that old F.R = new F
393.         if size of F = 2M x N, I on bottom, this returns its inverse in that range.
394.         basis = list of I-rows at end, index tells which column each corresponds to. """
395.         def gauss(F, M, N):
396.                 rno = len(F)
397.                 cno = len(F)
398.
399.                 index = [None]*M # eventually, gives column index of each basic variable
400.                 F = F.clone()
401.                 # now find an affine basis
402.                 m = 0; n = 0
403.                 others = range(M)
404.                 basis = []
405.                 while m < M and n < N:
406.                         # work with the matrix [m, M) x [n, N)
407.                         # locate first non-zero item in row[m]
408.                         row = F[m]
409.                         for j in range(n, N):
410.                                 if row[j] != 0:
411.                                         break
412.                         if j == N:
413.                                 m += 1
414.                                 break
415.                         else:
416.                                 # swap columns n, j in F and R
417.                                 # make 0 the rest of row[m]
418.                                 # and then divide column[n] by p = F[m][n]
419.                                 p = F[m][j]
420.                                 index[m] = len(basis)
421.                                 basis.append(m)
422.                                 others.remove(m)
423.                                 for i in range(rno):
424.                                         tmp = F[i][j]
425.                                         F[i][j] = F[i][n]
426.                                         F[i][n] = tmp/p
427.                                 # subtract column[a][m]*column[n] from column[a], all a != n
428.                                 for a in range(n):
429.                                         p = F[m][a]
430.                                         for b in range(rno):
431.                                                 F[b][a] -= p*F[b][n]
432.                                 for a in range(n+1, cno):
433.                                         p = F[m][a]
434.                                         for b in range(rno):
435.                                                 F[b][a] -= p*F[b][n]
436.                                 m += 1
437.                                 n += 1
438.                 return [F, index, basis, others]
439.         gauss=staticmethod(gauss)
440.
441. # -----------------------------------------------------------------------------------