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[0]
  32.                         self.basis = args[1]
  33.                         self.others = args[2]
  34.                         self.index = args[3]
  35.                         self.dim = None
  36.                         self.rank = len(F[0]) - 1
  37.                 else:
  38.                         N = len(F[0])-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
  92.                 # updates R also
  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[0])-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[0]) - 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[0])):
  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[0])):
  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[0]) - 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[0])-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[0]) - 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[0])
  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. # -----------------------------------------------------------------------------------

Raw Paste


Login or Register to edit or fork this paste. It's free.