(* Vertex Enumeration for Convex Polytopes and Arrangements Version 0.41 Beta November 14, 1992 Copyright (c) 1991,1992 by Komei Fukuda and Ichiro Mizukoshi This package contains Mathematica implementations of Avis-Fukuda algorithms for enumerating all vertices of a convex polyhedra given by a system of linear inequalities, and all vertices (points) of an arrangement given similarly. This package is copyright 1991 by Komei Fukuda and Ichiro Mizukoshi. This package may be copied in its entirety for nonprofit purposes only. Sale, other than for the direct cost of the media, is prohibited. This copyright notice must accompany all copies. The authors make no representations, express or implied, with respond to this documentation, of the software it describes and contains, including without limitations, any implied warranties of mechantability or fitness for a particular purpose, all of which are expressly disclaimed. The authors shall in no event be liable for any indirect, incidental, or consequential damages. This beta release is designed to run under Version 1.2 & 2.0 of Mathematica. Any comments, bug reports, or requests to get on the VertexEnumeration mailing list should be forwarded to: Komei Fukuda or Ichiro Mizukoshi Graduate School of Systems Management University of Tsukuba, Tokyo 3-29-1 Otsuka, Bunkyo-ku Tokyo 112, Japan +81-3-3942-6876 fukuda@gssm.otsuka.tsukuba.ac.jp or mizukosi@gssm.otsuka.tsukuba.ac.jp *) BeginPackage["VertexEnum`"] Unprotect[VE] Unprotect[PE] Unprotect[VertexEnumeration] Unprotect[PointEnumeration] Unprotect[CrissCrossSolve] Unprotect[BlandSolve] Unprotect[Faces] Unprotect[Polyhedron3D] Unprotect[SubsetQ] Unprotect[AdjacentQ] Unprotect[EdgesOfPolyhedron] Unprotect[MaximalQ] Unprotect[MaximalSets] Unprotect[LinesOfArrangement] VE::usage="VE[m,b] is a short form of VertexEnumeration[m,b]."; PE::usage="PE[m,b] is a short form of PointEnumeration[m,b]."; VertexEnumeration::usage="VertexEnumeration[m,b] gives the list of all vertices of the polyhedron represented as the solution set of the linear inequality system m.x <= b, x>=0. If you want to get results on progress MonitoringFile -> \"stdout\" "; PointEnumeration::usage="PointEnumeration[m,b] gives the list of all vertices of the arrangement of hyperplanes associated with the linear inequality system m.x <= b, x>=0."; CrissCrossSolve::usage="CrissCrossSolve[c,m,b] solves the linear program maximize c.x subject to m.x<=b and x>=0 by the Criss-Cross method. The output is {x*, v*, y*} where x* is an optimal solution, v* is an optimal value, and y* is a dual optimal solution. It warns if it is infeasible or dualinfeasible"; BlandSolve::usage="BlandSolve[c,m,b] solves the linear program maximize c.x subject to m.x<=b and x>=0 by the Simplex method with Bland's rule, provided b is nonnegative. The output is {x*, v*, y*} where x* is an optimal solution, v* is an optimal value, and y* is a dual optimal solution. It warns if it is unbounded."; Faces::usage="Faces[list] returns a list of faces for the zero-variable list. The kth face corresponds to kth inequality in x>=0, m.x<=b, for VertexEnumeration[m,b]."; Polyhedron3D::usage="Polyhedron3D[ver_List,zerova_List] returns 3D graphic primitives for drawing a 3D polytope. The two arguments correspond to the first two list given by VertexEnumeration[m,b]."; SubsetQ::usage="SubsetQ[s,t] returns true if s is a subset of t." AdjacentQ::usage="AdjacentQ[i_Integer,j_Integer,zerova_List] returns True if and only if the zeroset zerova[[i]] and the zeroset zerova[[j]] correspond to adjacent vertices (the intersection of zerova[[i]] and zerova[[j]] is contained in no other zerosets)." EdgesOfPolyhedron::usage="EdgesOfPolyhedron[zerova_List] returns the set of adjacent pairs of vertices. See AdjacentQ." MaximalQ::usage="MaximalQ[i_Integer,l_List] returns True if and only if l[[i]] is a maximal set." MaximalSets::usage="MaximalSets[l_List] returns the sublist of maximal sets. See, MaximalQ." LinesOfArrangement::usage="LinesOfArrangement[zerova_List] returns the list of lines as [1] sets of colinear vertices, and as [2] zerosets." Options[VertexEnumeration] = {SearchTree -> False, ZeroVariables -> True, MonitoringFile -> "/dev/null"} Options[PointEnumeration] = {SearchTree -> False, ZeroVariables -> True, MonitoringFile -> "/dev/null"} Begin["`Private`"] VertexEnumeration::LpInfeasible = "linear program is infeasible." VertexEnumeration::LpDualInfeasible = "linear program is dual infeasible." VertexEnumeration::LpUnbounded = "linear program is unbounded." VertexEnumeration::Infeasible = "linear inequality system is infeasible." (* comments for private routine *) Pivot::usage="Pivot[mat,r,s,b,n] performs a pivot operation on mat with row:r and column:s and updates basis b and nonbasis n. Last 2 arguments are optional."; GetPivotSequence::usage="GetPivotSequence[startbv,startnbv,goalbv,goalnbv] gives a set of pivots (variable pairs) which transfers the start to the goal."; SortDict::usage="SortDict[mat,b,n] returns the new {mat,b,n} sorted by b and n reordered in ascending order."; Dic2Vertex::usage="Dic2Vertex[value of bv, name of bv, name of nbv, ver] returns the vertex associated with ver."; MakeSequence::usage="MakeSequence[dic_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,bases_?MatrixQ] returns pivot sequence."; GetChangeBaseVariableNames::usage="GetChangeBaseVariableNames[startbv,startnbv,goalbv,goalnbv] gives a set of pivots (variable pairs) which transfers the start to the goal."; ExecPivotSequence::usage="ExecPivotSequence[dic,bv,nbv,sequence] performs a sequence of pivots to {dic,bv,nbv} according to the given sequence." NonEmptyFaces::usage="NonEmptyFaces[list] returns a list of nonempty faces."; OptimalBasesEnumeration::usage="OptimalBasesEnumeration[mat,b] gives the set of variable pairs obtaining an optimal basis."; LexMin::usage="LexMin[mat,bv,nbv,m,n] returns True if the dictionary {mat,bv,nbv} is lexicographically minimum for representing the current solution."; CrissCrossSelect::usage="CrissCrossSelect[mat,bv,nbv,m,n] gives the pivot position selected by the Criss-Cross rule."; CrissCrossSelectQ::usage="CrissCrossSelectQ[mat,bv,nbv,i,j] checks whether {i,j} is the pivot position selected by the Criss-Cross rule."; CrissCrossFilter::usage="CrissCrossFilter[mat,bv,nbv,i,j] is filter (necesarily True) for {i,j}==CrissCrossSelect[mat,bv,nbv]."; CrissCross::usage="CrissCross[mat,bv,nbv,m,n] solves the linear program given by a dictionary {mat,bv,nbv} and outputs a terminal dictionary."; NegSmallest::usage="Ratio Function used by the Simplex Method"; FirstElement::usage="FirstElement[function,list1[,list2]] gives a position on list1 (ordered by list2),that is True at first time."; NonNegativeVectorQ::usage="NonNegativeVectorQ[vec] returns True if and only if vec is a nonnegative vector."; BlandSelect::usage="BlandSelect[mat,bv,nbv] gives a pivot position selected by Bland's rule."; BlandSelectQ::usage="BlandSelectQ[mat,bv,nbv,i,j] returns True if and only if {i,j}== BlandSelect[mat,bv,nbv]."; BlandFilter::usage="BlandFilter[mat,bv,nbv,i,j,m,n] is filter (necesarily True) for {i,j}== BlandSelect[mat,bv,nbv]."; Bland::usage="Bland[mat,bv,nbv,m,n] solves the linear program given by a feasible dictionary {mat,bv,nbv} and outputs a terminal dictionary."; DegenerateQ::usage="DegenerateQ[vec] gives True if vec contains 0. It is used to recognize degenerate dictionary."; InfeasibleQ::usage="InfeasbleQ[mat] gives True if mat is infeasible."; DualInfeasibleQ::usage="DualInfeasbleQ[mat] gives True if mat is dual infeasible."; Zerovar::usage="Zerovar[Value of bv, Name of bv,Name of nbv] shows the variable names with value 0." BSearch::usage="BSearch[mat,b,n] returns the set of all basic feasible solutions for the input mat, b, n, {value of bv, name of bv, name of nbv}."; CSearch::usage="CSearch[mat,b,n] returns the set of all basic solutions for the input mat, b, n, {value of bv, name of bv, name of nbv}."; Options[CrissCross] = { NegativeTest-> Function[Negative[Chop[N[Simplify[#]]]]]} Options[Bland] = { NegativeTest-> Function[Negative[Chop[N[Simplify[#]]]]]} Options[ZeroVar] = { NegativeTest-> Function[Negative[Chop[N[Simplify[#]]]]]} Options[BSearch] = {SearchTree -> False, NegativeTest-> Function[Negative[Chop[N[Simplify[#]]]]], MonitoringFile->"/dev/null"} Options[CSearch] = {SearchTree -> False, NegativeTest-> Function[Negative[Chop[N[Simplify[#]]]]], MonitoringFile->"/dev/null"} VE[mat_?MatrixQ,vec_?VectorQ, opts___Rule]:=VertexEnumeration[mat,vec,opts] PE[mat_?MatrixQ,vec_?VectorQ, opts___Rule]:=PointEnumeration[mat,vec,opts] Bout2Dlist[bout_List,dlist0_List]:= Block[{i,j,dlist,boutlength,dlist0rule}, dlist = Map[ Function[Father2Edge[Dlist2Father[#[[2]]]]], bout]; boutlength = Map[Function[Length[#[[1]]]],bout]; boutlength = Prepend[Table[Sum[boutlength[[i]],{i,1,j}],{j,Length[boutlength]-1}],0]; dlist += boutlength; dlist0rule = Table[ Rule[i, boutlength[[i]] + 1], {i,Length[boutlength]}]; AppendTo[dlist,dlist0 /. dlist0rule]; {Flatten[dlist,1],boutlength+1} ] VertexEnumeration[mat_?MatrixQ,vec_?VectorQ,opts___Rule]:= Block[{dic,i,j,m,n,bv,nbv,flattenbout,bout,obe,output,dlist,opt1,opt2}, opt1 = SearchTree /. {opts} /. Options[VertexEnumeration]; opt2 = ZeroVariables /. {opts} /. Options[VertexEnumeration]; {m,n} = Dimensions[mat]; dic = Transpose[Append[Transpose[-mat],vec]]; AppendTo[dic,Append[Table[0,{n}],0]]; (* get dic, bv, nbv for Bland. Warnig nbv has been chanded.*) {dic,bv,nbv} = CrissCross[dic,Range[n+1,m+n],Range[n],m+1,n+1,opts]; (* feasiblity check *) If[InfeasibleQ[dic], Message[VertexEnumeration::Infeasible];Return[{}]]; (* shape up dic for Bland *) dic = Drop[dic,-1]; If[DegenerateQ[Last[Transpose[dic]]], (* There is degeneracy *) (* obe looks like {{{dic,bv,nbv},{dic,bv,nbv}}, dlist0} *) obe= OptimalBasesEnumeration[Table[-dic[[i,j]],{i,m},{j,n}],Table[dic[[i,n+1]],{i,m}],bv,nbv]; (* bout looks like {{{{Value of bv,bv,nbv},{Value of bv,bv,nbv},...},dlist},{{{Value of bv,bv,nbv},...},dlist},...} *) bout = Map[ Function[BSearch[#[[1]],#[[2]],#[[3]],opts]],obe[[1]]]; (* flattenbout looks like {{Value of bv,bv,nbv},{Value of bv,bv,nbv},....} but no dupilicat startpoint*) flattenbout = Flatten[Map[Function[#[[1]]],bout],1]; output = {Map[Function[Dic2Vertex[ #[[1]],#[[2]],#[[3]],Range[n]]],flattenbout]}; If[opt2, AppendTo[output, Map[ Function[ZeroVar[#[[1]],#[[2]],#[[3]],opts]], flattenbout]]]; If[opt1, dlist = Bout2Dlist[bout,obe[[2]]]; AppendTo[output,dlist[[1]]]; AppendTo[output,dlist[[2]]]], (* There is not degeneracy *) (* dic and mat must be same shape *) AppendTo[dic,Append[Table[-1,{n}],0]]; bout = BSearch[dic,bv,nbv,opts]; output = {Map[ Function[Dic2Vertex[ #[[1]],#[[2]],#[[3]],Range[n]]],bout[[1]]]}; If[opt2, AppendTo[output, Map[ Function[ZeroVar[#[[1]],#[[2]],#[[3]],opts]], bout[[1]]]]]; If[opt1, AppendTo[output, Father2Edge[Dlist2Father[bout[[2]]]]]] ]; output ] PointEnumeration[mat_?MatrixQ,vec_?VectorQ,opts___Rule]:= Block[{dic,i,j,m,n,bv,nbv,tmp,flattenbout,bout,obe,output,opt1,opt2}, opt1 = SearchTree /. {opts} /. Options[PointEnumeration]; opt2 = ZeroVariables /. {opts} /. Options[PointEnumeration]; {m,n} = Dimensions[mat]; (* dic = Table[ Which[ SimpleNegative[vec[[i]]], Join[mat[[i]],vec[[i]]], True, Join[-mat[[i]],-vec[[i]]]], {i,m} ] *) dic = Transpose[Append[Transpose[-mat],vec]]; AppendTo[dic,Append[Table[-1,{n}],0]]; bv = Range[n+1,m+n]; nbv = Range[n]; If[DegenerateQ[Last[Transpose[dic]]], (* There is degeneracy *) (* obe looks like {{{dic,bv,nbv},{dic,bv,nbv}}, dlist0} *) obe= OptimalBasesEnumeration[Table[-dic[[i,j]],{i,m},{j,n}],Table[dic[[i,n+1]],{i,m}],bv,nbv]; (* bout looks like {{{{Value of bv,bv,nbv},{Value of bv,bv,nbv},...},dlist},{{{Value of bv,bv,nbv},...},dlist},...} *) bout = Map[ Function[CSearch[#[[1]],#[[2]],#[[3]],opts]],obe[[1]]]; (* flattenbout looks like {{Value of bv,bv,nbv},{Value of bv,bv,nbv},....} but no dupilicat startpoint*) flattenbout = Flatten[Map[Function[#[[1]]],bout],1]; output = {Map[Function[Dic2Vertex[ #[[1]],#[[2]],#[[3]],Range[n]]],flattenbout]}; If[opt2, AppendTo[output, Map[ Function[ZeroVar[#[[1]],#[[2]],#[[3]],opts]], flattenbout]]]; If[opt1, dlist = Bout2Dlist[bout,obe[[2]]]; AppendTo[output,dlist[[1]]]; AppendTo[output,dlist[[2]]]], (* There is not degeneracy *) (* dic and mat must be same shape *) tmp = CSearch[dic,bv,nbv,opts]; output = {Map[ Function[Dic2Vertex[ #[[1]],#[[2]],#[[3]],Range[n]]],tmp[[1]]]}; If[opt2, AppendTo[output, Map[ Function[ZeroVar[#[[1]],#[[2]],#[[3]],opts]], tmp[[1]]]]]; If[opt1, AppendTo[output, Father2Edge[Dlist2Father[tmp[[2]]]]]] ]; output ] Faces[zerova_List]:= Map[Function[x, If[Position[zerova,x]=={},{}, Transpose[Position[zerova,x]][[1]] ] ], Range[Max[zerova]] ] NonEmptyFaces[zerova_List]:= Map[Function[Transpose[Position[zerova,#]][[1]]],Union[Flatten[zerova]]] MakeLink[ver_List,zerova_List]:= Block[{candidates, cycle ,candidate2cycle}, candidates = Drop[ver,1]; cycle = { First[ver] }; (* when one is list, comparing needs other is list *) While[{} != (candidate2cycle = Select[candidates,AdjacentQ[Last[cycle],#,zerova] &]), AppendTo[cycle,candidate2cycle[[1]]]; candidates = Complement[candidates,candidate2cycle[[{1}]]]; ]; While[{} != (candidate2cycle = Select[candidates,AdjacentQ[First[cycle],#,zerova] &]), PrependTo[cycle,candidate2cycle[[1]]]; candidates = Complement[candidates,candidate2cycle[[{1}]]]; ]; cycle ] Polyhedron3D[ver_List,zerova_List]:= Polygon /@ (ver[[#]]&) /@ ( MakeLink[#,zerova] & ) /@ NonEmptyFaces[zerova] ZeroVar[ValueOfBv_?VectorQ,bv_?VectorQ,nbv_?VectorQ,opts___Rule]:= Block[{m,out={}}, m=Length[ValueOfBv]; SimpleNegative := NegativeTest /. {opts} /. Options[ZeroVar]; out=Map[Function[SimpleSameQ[#,0]],ValueOfBv]; Join[nbv,Map[Function[bv[[#[[1]]]]],Position[out,True]]] ] ForAllQ[func_,list1_List]:= Apply[And, Map[func, list1]] Father2Edge[fathers_?VectorQ]:= Map[{#,fathers[[#]]}&,Range[2,Length[fathers]]] Dlist2Father[dlist_?VectorQ]:= Block[{po,dist,father,fathers={1}}, (* first element in dist must be 0 and Length[dlist] => 2*) For[po = 2, po <= Length[dlist], ++po, dist = dlist[[po]]; father = po - 1; --dist; (* At first step a element in dist must be 1 *) While[dist > 0, father = fathers[[father]]; --dist ]; AppendTo[fathers,father]; ]; fathers ] Dic2Vertex[valueofbv_?VectorQ,bv_?VectorQ,nbv_?VectorQ,ver_?VectorQ]:= Map[Function[If[ MemberQ[nbv,#], 0,valueofbv[[Position[bv,#][[1,1]]]]]], ver] LexMin[dic_?MatrixQ,bv_?VectorQ,nbv_?VectorQ]:= Block[{m,n}, {m,n}=Dimensions[dic]; LexMin[dic,bv,nbv,m,n] ] LexMin[dic_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,m_Integer,n_Integer]:= Block[{i,j}, For[i = 1, i < m , ++i, If[SimpleSameQ[0,dic[[i,n]]], For[j = 1, j < n, ++j, If[ Not[SimpleSameQ[0,dic[[i,j]]]] && bv[[i]] > nbv[[j]], Return[False]] ] ] ]; True ] CSearch[dict_?MatrixQ,opts___Rule]:= Block[{m,n,bv,nbv}, {m,n}=Dimensions[dict]; bv = Range[m-1]; nbv = Range[m,m+n-2]; CSearch[dict,bv,nbv,opts] ] CSearch[origdict_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,opts___Rule]:= Block[{ bv=origbv,nbv=orignbv,dict=origdict,m,n,i=1,j=1,opt1,MonitoringFile, result,dist=1,distlist={0}}, opt1 = SearchTree /. {opts} /. Options[CSearch]; MonitoringFile = MonitoringFile /. {opts} /. Options[CSearch]; SimpleNegative := NegativeTest /. {opts} /. Options[CSearch]; {m,n}=Dimensions[dict]; If[opt1, result ={{Drop[Transpose[origdict][[n]],-1],origbv,orignbv}}; Write[MonitoringFile,Last[result][[3]]], If[LexMin[origdict,origbv,orignbv,m,n], result ={{Drop[Transpose[origdict][[n]],-1],origbv,orignbv}}; Write[MonitoringFile,Last[result][[3]]], result = {}]]; While[ ((i < m )|| (origbv != bv)), While[ ((i < m ) && Not[CrissCrossSelectQ[dict,bv,nbv,i,j]]), ++j; If[ j >= n, j = 1; ++i]; ]; If[ i < m, {dict,bv,nbv} = Pivot[dict,bv,nbv,i,j,m,n]; If[opt1, AppendTo[result,{Drop[Transpose[dict][[n]],-1],bv,nbv}]; Write[MonitoringFile,Last[result][[3]]]; AppendTo[distlist,dist]; dist = 1, If[LexMin[dict,bv,nbv,m,n], AppendTo[result,{Drop[Transpose[dict][[n]],-1],bv,nbv}]; Write[MonitoringFile,Last[result][[3]]]; AppendTo[distlist,dist]; dist =1; ]; ]; {i,j}={1,1}, {i,j}=CrissCrossSelect[dict,bv,nbv,m,n]; If[ (i < m && j < n), {dict,bv,nbv}=Pivot[dict,bv,nbv,i,j,m,n]; ++j; ++dist; If[ j >= n, j = 1; ++i] ]; ]; ]; {result,distlist} ] BSearch[dict_?MatrixQ,opts___Rule]:= Block[{m,n,bv,nbv}, {m,n}=Dimensions[dict]; bv = Range[m-1]; nbv = Range[m, m+n-2]; BSearch[dict,bv,nbv,opts] ] BSearch[origdict_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,opts___Rule]:= Block[{ bv=origbv,nbv=orignbv,dict=origdict,m,n,i=1,j=1,opt1,MonitoringFile, result,dist=1,distlist = {0},bvposi,nbvposi}, opt1 = SearchTree /. {opts} /. Options[BSearch]; MonitoringFile = MonitoringFile /. {opts} /. Options[BSearch]; SimpleNegative := NegativeTest /. {opts} /. Options[BSearch]; {m,n}=Dimensions[dict]; bvposi = Flatten[Map[Position[bv,#] & ,Sort[bv]]]; nbvposi = Flatten[Map[Position[nbv,#] & ,Sort[nbv]]]; If[opt1, result ={{Drop[Transpose[origdict][[n]],-1],origbv,orignbv}}; Write[MonitoringFile,Last[result][[3]]], If[LexMin[origdict,origbv,orignbv,m,n], result ={{Drop[Transpose[origdict][[n]],-1],origbv,orignbv}}; Write[MonitoringFile,Last[result][[3]]], result = {}]]; While[ ((i < m) || (origbv != bv)), While[ (i < m && Not[BlandSelectQ[dict,bv,nbv,i,j,m,n,bvposi]]), ++j; If[ j >= n, j = 1; ++i] ]; If[ i < m, {dict,bv,nbv} = Pivot[dict,bv,nbv,i,j,m,n]; bvposi = Flatten[Map[Position[bv,#] & ,Sort[bv]]]; nbvposi = Flatten[Map[Position[nbv,#] & ,Sort[nbv]]]; (* Print["."]; *) If[opt1, AppendTo[result,{Drop[Transpose[dict][[n]],-1],bv,nbv}]; Write[MonitoringFile,Last[result][[3]]]; AppendTo[distlist,dist]; dist = 1, If[LexMin[dict,bv,nbv,m,n], AppendTo[result,{Drop[Transpose[dict][[n]],-1],bv,nbv}]; Write[MonitoringFile,Last[result][[3]]]; AppendTo[distlist,dist]; dist = 1, ]; ]; {i,j}={1,1}, {i,j}=BlandSelect[dict,bv,nbv,bvposi,nbvposi]; If[ (i < m && j < n), {dict,bv,nbv}=Pivot[dict,bv,nbv,i,j,m,n]; bvposi = Flatten[Map[Position[bv,#] & ,Sort[bv]]]; nbvposi = Flatten[Map[Position[nbv,#] & ,Sort[nbv]]]; ++j; ++dist; If[ j >= n, j = 1; ++i] ]; ]; ]; {result,distlist} ] BlandFilter[dict_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,i_Integer,j_Integer]:= Block[{m,n}, {m,n}=Dimensions[dict]; BlandFilter[dict,bv,nbv,i,j,m,n,Flatten[Map[Position[bv,#] & ,Sort[bv]]]] ] BlandFilter[dict_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,i_Integer,j_Integer,m_Integer,n_Integer]:= BlandFilter[dict,bv,nbv,i,j,m,n,Flatten[Map[Position[bv,#] & ,Sort[bv]]]] BlandFilter[dict_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,i_Integer,j_Integer,m_Integer,n_Integer,bvposi_?VectorQ]:= ( SimpleNegative[dict[[i,j]]] && SimpleNegative[dict[[m,j]]] && SimpleSameQ[-dict[[i,n]]/dict[[i,j]], NegSmallest[Transpose[dict][[j]],bvposi,Last[Transpose[dict]]][[2]]] ) BlandSelectQ[origdict_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,i_Integer,j_Integer,m_Integer,n_Integer]:= BlandSelectQ[origdict,origbv,orignbv,i,j,m,n,Flatten[Map[Position[origbv,#] & ,Sort[origbv]]]] BlandSelectQ[origdict_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,i_Integer,j_Integer,m_Integer,n_Integer,bvposi_?VectorQ]:= Block[{dict,bv,nbv}, BlandFilter[origdict,origbv,orignbv,i,j,m,n,bvposi] && ( {dict,bv,nbv}=Pivot[origdict,origbv,orignbv,i,j,m,n]; SameQ[{i,j},BlandSelect[dict,bv,nbv]] ) ] NonNegativeVectorQ[vec_?VectorQ]:= ForAllQ[Function[Not[SimpleNegative[#]]],vec] BlandSolve[cost_?VectorQ,mat_?MatrixQ,b_?NonNegativeVectorQ,opts___Rule]:= Block[{dic,m,n,ccdic,ccbv,ccnbv}, dic = Transpose[Append[Transpose[Append[-mat,cost]],Append[b,0]]]; {m,n}=Dimensions[dic]; {ccdic, ccbv, ccnbv}= Bland[dic,opts]; If[DualInfeasibleQ[ccdic], Message[VertexEnumeration::LpUnbounded];Return[{}]]; {Dic2Vertex[ Drop[Transpose[ccdic][[n]],-1], ccbv, ccnbv,Range[n-1]],ccdic[[m,n]], Dic2Vertex[ Drop[-ccdic[[m]],-1],ccnbv, ccbv, Range[n,n+m-2]]} ] Bland[dict_?MatrixQ,opts___Rule]:= Block[{m,n,bv,nbv}, {m,n}=Dimensions[dict]; bv = Range[n,n+m-2]; nbv = Range[n-1]; Bland[dict,bv,nbv,m,n,opts] ] Bland[origdict_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,opts___Rule]:= Block[{m,n}, {m,n}=Dimensions[origdict]; Bland[origdict,origbv,orignbv,m,n,opts] ] Bland[origdict_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,m_Integer,n_Integer,opts___Rule]:= Block[{bv=origbv,nbv=orignbv,dict=origdict,i,j}, SimpleNegative := NegativeTest /. {opts} /. Options[Bland]; For[{i,j}=BlandSelect[dict,bv,nbv], (i < m && j < n), {i,j}=BlandSelect[dict,bv,nbv], {dict,bv,nbv}=Pivot[dict,bv,nbv,i,j,m,n] ]; {dict,bv,nbv} ] BlandSelect[dict_?MatrixQ,bv_?VectorQ,nbv_?VectorQ]:= Block[{bvposi,nbvposi}, bvposi = Flatten[Map[Position[bv,#] & , Sort[bv]]]; nbvposi = Flatten[Map[Position[nbv,#] & , Sort[nbv]]]; BlandSelect[dict,bv,nbv,bvposi,nbvposi] ] BlandSelect[dict_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,bvposi_?VectorQ,nbvposi_?VectorQ]:= Block[{i,j}, j = FirstElement[SimplePositive,Last[dict],nbvposi]; i = NegSmallest[ Transpose[dict][[j]],bvposi, Last[Transpose[dict]] ][[1]]; {i,j} ] CrissCrossFilter[dict_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,i_Integer,j_Integer]:= Block[{m,n}, {m,n}=Dimensions[dict]; CrissCrossFilter[dict,bv,nbv,i,j,m,n] ] CrissCrossFilter[dict_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,i_Integer,j_Integer,m_Integer,n_Integer]:= (( SimplePositive[dict[[i,n]]] && SimplePositive[dict[[i,j]]] && ForAllQ[Function[Not[SimpleNegative[dict[i,#]]] && nbv[#] < bv[i]],Range[n-1]] ) || ( SimpleNegative[dict[[m,j]]] && SimpleNegative[dict[[i,j]]] && ForAllQ[Function[Not[SimplePositive[dict[#,j]]] && bv[#] < nbv[j]],Range[m-1]] )) CrissCrossSelectQ[origdict_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,i_Integer,j_Integer]:= Block[{dict,bv,nbv,m,n}, {m,n}=Dimensions[origdict]; CrissCrossFilter[origdict,origbv,orignbv,i,j,m,n] && ( {dict,bv,nbv}=Pivot[origdict,origbv,orignbv,i,j,m,n]; SameQ[{i,j},CrissCrossSelect[dict,bv,nbv,m,n]] ) ] CrissCrossSolve[cost_?VectorQ,mat_?MatrixQ,b_?VectorQ,opts___Rule]:= Block[{dic,m,n,ccdic,ccbv,ccnbv}, dic = Transpose[Append[Transpose[Append[-mat,cost]],Append[b,0]]]; {m,n}=Dimensions[dic]; {ccdic, ccbv, ccnbv}=CrissCross[dic,opts]; If[InfeasibleQ[ccdic], Message[VertexEnumeration::LpInfeasible];Return[{}]]; If[DualInfeasibleQ[ccdic], Message[VertexEnumeration::LpDualInfeasible];Return[{}]]; {Dic2Vertex[ Drop[Transpose[ccdic][[n]],-1], ccbv, ccnbv,Range[n-1]],ccdic[[m,n]], Dic2Vertex[ Drop[-ccdic[[m]],-1],ccnbv, ccbv, Range[n,n+m-2]]} ] CrissCross[dict_?MatrixQ,opts___Rule]:= Block[{m,n,bv,nbv}, {m,n}=Dimensions[dict]; bv = Range[n,m+n-2]; nbv = Range[n-1]; CrissCross[dict,bv,nbv,m,n,opts] ] CrissCross[origdict_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,opts___Rule]:= Block[{m,n}, {m,n}=Dimensions[origdict]; CrissCross[origdict,origbv,orignbv,m,n,opts] ] CrissCross[origdict_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,m_Integer,n_Integer,opts___Rule]:= Block[{bv=origbv,nbv=orignbv,dict=origdict,i,j}, SimpleNegative := NegativeTest /. {opts} /. Options[CrissCross]; For[{i,j}=CrissCrossSelect[dict,bv,nbv,m,n],(i < m && j < n), {i,j}=CrissCrossSelect[dict,bv,nbv,m,n], {dict,bv,nbv}=Pivot[dict,bv,nbv,i,j,m,n] ]; {dict,bv,nbv} ] CrissCrossSelect[dict_?MatrixQ,bv_?VectorQ,nbv_?VectorQ]:= Block[{m,n}, {m,n}=Dimensions[dict]; CrissCrossSelect[dict,bv,nbv,m,n] ] CrissCrossSelect[dict_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,m_Integer,n_Integer]:= Block[{i,j,bvposi,nbvposi}, bvposi = Flatten[Map[ Position[bv,#] & ,Sort[bv]]]; nbvposi = Flatten[Map[ Position[nbv,#] & ,Sort[nbv]]]; i = FirstElement[SimpleNegative,Last[Transpose[dict]],bvposi]; j = FirstElement[SimplePositive,Last[dict],nbvposi]; Which[ i == m, i = FirstElement[SimpleNegative,Transpose[dict][[j]],bvposi], j == n, j = FirstElement[SimplePositive,dict[[i]],nbvposi], nbv[[j]] > bv[[i]], j = FirstElement[SimplePositive,dict[[i]],nbvposi], True, i = FirstElement[SimpleNegative,Transpose[dict][[j]],bvposi]; ]; {i,j} ] FirstElement[func_,list1_List,list2_List]:= Block[{i}, For[i = 1,i <= Length[list2],++i, If[func[list1[[list2[[i]]]]],Return[list2[[i]]]]]; i ] FirstElement[func_,list1_List]:= Block[{i}, For[i = 1, i <= Length[list1], ++i, If[func[list1[[i]]],Return[i]]]; i ] NegSmallest[list1_List,list2_List,list3_List]:= Block[{i,posi,tmpval,val=Infinity}, For[i = 1,i <= Length[list2],++i, If[ ( SimpleNegative[list1[[list2[[i]]]]] && SimpleGreater[val, tmpval = -(list3[[list2[[i]]]]/list1[[list2[[i]]]])] ), posi = i; val = tmpval ]; ]; (* val = Infinity -> nothing happens. Never happens *) If[SameQ[Infinity,val], {i,val}, {list2[[posi]],val}] ] Pivot[d_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,i_Integer,j_Integer,m_Integer,n_Integer]:= Block[{tmpbv,tmpnbv}, {tmpbv,tmpnbv} = {bv,nbv}; {tmpbv[[i]],tmpnbv[[j]]}={tmpnbv[[j]],tmpbv[[i]]}; {Pivot[d,i,j,m,n],tmpbv,tmpnbv} ] Pivot[d_?MatrixQ,bv_?VectorQ,nbv_?VectorQ,i_Integer,j_Integer]:= Block[{tmpbv,tmpnbv}, {tmpbv,tmpnbv} = {bv,nbv}; {tmpbv[[i]],tmpnbv[[j]]}={tmpnbv[[j]],tmpbv[[i]]}; {Pivot[d,i,j],tmpbv,tmpnbv} ] Pivot[d_?MatrixQ,r_Integer,s_Integer]:= Block[{m,n}, {m,n}=Dimensions[d]; Pivot[d,r,s,m,n] ] Pivot[d_?MatrixQ,r_Integer,s_Integer,m_Integer,n_Integer]:= Block[{i,j}, Simplify[ Table[ Which[ (i == r) && (j == s), 1 / d[[i,j]], (i == r) && (j != s), -(d[[i,j]]/d[[r,s]]), (i != r) && (j == s), d[[i,j]]/d[[r,s]], True, d[[i,j]] - ((d[[r,j]] * d[[i,s]]) / d[[r,s]]) ],{i,m},{j,n} ] ] ] GetChangeBaseVariableNames[startbv_?VectorQ,startnbv_?VectorQ,goalbv_?VectorQ,goalnbv_?VectorQ]:= {Intersection[startbv,goalnbv],Intersection[startnbv,goalbv]} MakeSequence[origdic_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,bases_?MatrixQ]:= Block[{dic = origdic, bv = origbv, nbv = orignbv, bvposi, nbvposi, tmpsequence}, (* name to position *) bvposi = Flatten[Map[Position[bv,#] & ,bases[[1]]]]; nbvposi = Flatten[Map[Position[nbv,#] & ,bases[[2]]]]; (* make all possible sequences *) tmpsequence = Map[Transpose[{bvposi,#}] &, Permutations[nbvposi]]; (* choose one *) If[{} == tmpsequence,tmpsequence,tmpsequence[[FirstElement[# &,Map[Apply[And,Map[Not[SimpleSameQ[0,origdic[[#[[1]],#[[2]]]]]] &,#]] &,tmpsequence]]]]] ] GetPivotSequence[startbv_?VectorQ,startnbv_?VectorQ,goalbv_?VectorQ,goalnbv_?VectorQ]:= Transpose[{Intersection[startbv,goalnbv],Intersection[startnbv,goalbv]}] ExecPivotSequence[origdic_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ,sequence_?MatrixQ]:= Block[{dic = origdic, bv = origbv, nbv = orignbv, i}, Do[{dic,bv,nbv} = Pivot[dic,bv,nbv,sequence[[i,1]],sequence[[i,2]]],{i,Length[sequence]}]; {dic,bv,nbv} ] OptimalBasesEnumeration[mat_?MatrixQ,vec_?VectorQ,bv_?VectorQ,nbv_?VectorQ]:= Block[{dic,degenerateBasicVariables,i,m,n,tmp,output}, {m,n} = Dimensions[mat]; degenerateBasicVariables = Flatten[Position[vec,0]]; dic = mat[[degenerateBasicVariables]]; AppendTo[dic,Table[1,{n}]]; dic = Append[Transpose[dic],Append[Table[-1,{Length[degenerateBasicVariables]}],0]]; (* bv = Table[i,{i,n+1,m+n}]; nbv = Table[i,{i,n}]; *) (* this is dual, so bv and nbv are exchanged *) tmp = BSearch[dic,nbv,bv[[degenerateBasicVariables]],SearchTree->True,MonitoringFile->"/dev/null"]; dic = Transpose[Append[Transpose[-mat],vec]]; AppendTo[dic,Append[Table[-1,{n}],0]]; output = Map[ Function[MakeSequence[dic,bv,nbv,GetChangeBaseVariableNames[bv,nbv,#[[3]],#[[2]]]]], tmp[[1]]]; output = {Map[Function[ExecPivotSequence[dic,bv,nbv,#]], output]}; AppendTo[output, Father2Edge[Dlist2Father[tmp[[2]]]]]; output ] SubsetQ[s_List,t_List]:= Length[s]==Length[Intersection[s,t]]; AdjacentQ[i_Integer,j_Integer,l_List]:= Block[{subface}, subface=Intersection[l[[i]],l[[j]]]; Length[Select[l,SubsetQ[subface,#]&]]==2] EdgesOfPolyhedron[l_List]:= Block[{i,j,edges={}}, Do[ Do[If[AdjacentQ[i,j,l],AppendTo[edges,{i,j}]], {j,i+1,Length[l]} ],{i,Length[l]} ]; edges] MaximalQ[i_Integer,l_List]:= Block[{candidate}, candidate=l[[i]]; Length[Union[Select[l,SubsetQ[candidate,#]&]]]==1] MaximalSets[l_List]:= Block[{i,maximals={}}, Do[If[MaximalQ[i,l], AppendTo[maximals,l[[i]]] ], {i,Length[l]} ]; maximals] LinesOfArrangement[zerova_List]:= Block[{i,j,p=Length[zerova],zeros={},lines={}}, Do[Do[AppendTo[zeros, Intersection[zerova[[i]],zerova[[j]]]],{j,i+1,p} ],{i,p} ];zeros=MaximalSets[Union[zeros]]; Do[AppendTo[lines, Select[Range[p],SubsetQ[zeros[[i]],zerova[[#]]]&] ],{i,Length[zeros]} ]; {lines,zeros} ] SortDict[origdict_?MatrixQ,origbv_?VectorQ,orignbv_?VectorQ]:= Block[{bv=origbv,nbv=orignbv,dict=origdict}, dict = Append[ Map[ Function[dict[[Position[bv,#][[1,1]]]]], Sort[bv]], Last[dict]]; dict = Append[ Map[ Function[Transpose[dict][[Position[nbv,#][[1,1]]]]], Sort[nbv]], Last[Transpose[dict]]]; {Transpose[dict],Sort[bv],Sort[nbv]} ] DegenerateQ[vec_?VectorQ]:=Count[vec,0]!=0 InfeasibleRowQ[vec_?VectorQ]:= SimpleNegative[Last[vec]] && ForAllQ[Function[Not[SimplePositive[#]]],Drop[vec,-1]] InfeasibleQ[mat_?MatrixQ]:= Apply[Or,Map[InfeasibleRowQ,Drop[mat,-1]]] DualInfeasibleQ[mat_?MatrixQ]:= InfeasibleQ[Transpose[-mat]] SimplePositive:=Function[SimpleNegative[-#]] SimpleSameQ:= Function[( Not[SimpleNegative[#1 - #2]] && Not[SimplePositive[#1 - #2]] )] SimpleGreater:=Function[SimplePositive[#1 - #2]] End[(* "`Private`" *)] Protect[VE] Protect[PE] Protect[VertexEnumeration] Protect[PointEnumeration] Protect[CrissCrossSolve] Protect[BlandSolve] Protect[Faces] Protect[Polyhedron3D] Protect[SubsetQ] Protect[AdjacentQ] Protect[EdgesOfPolyhedron] Protect[MaximalQ] Protect[MaximalSets] Protect[LinesOfArrangement] EndPackage[]