-- $Id: intprog.cpkg,v 4.6 2004/11/25 13:58:01 cocoa Exp $ Package $contrib/intprog -- Integer Programming Alias IP := $contrib/intprog, Mat := $mat; Define About() PrintLn " KeyWords : Integer Programming, Test Set, Cost Author : A.M.Bigatti Version : CoCoA3.7 Date : 18 June 1999 " EndDefine; -- About ------[ Manual ]-------- Define Man() PrintLn " Suggested alias for this package: Alias IP := $contrib/intprog; SYNTAX IP.TestSet(A: MAT, Cost: LIST): TAGGED('IP') IP.TestSet(IP: TAGGED('IP'), Cost: LIST): TAGGED('IP') IP.Sol(A: MAT, B: LIST, Cost: LIST); IP.Sol(IP: TAGGED('IP'), B: LIST, Cost: LIST); DESCRIPTION Let A be an NxM matrix with integer entries, B a list with N integer entries, C a list with M rational entries defining a linear cost function IP.TestSet(A, Cost) computes the toric ideal associated to A and a GBasis wrt a special order. This GBasis is a Test-Set for the Integer Programming problem IP_c. IP.TestSet(IP, Cost) same as IP.TestSet(A, Cost) using the precomputed information in IP IP.Sol(A, B, Cost); finds a minimal solution to AX = B wrt to the cost function IP.Sol(IP, B, Cost); same as IP.Sol(A, B, Cost) using the precomputed information in IP >EXAMPLE< Use R ::= Q[xyz]; -- any ring -- input with NON-NEGATIVE entries A := Mat[ [3, 1, 11, 2, 3, 5, 3], [4, 5, 0, 1, 7, 4, 6], [5, 6, 1, 9, 2, 3, 3]]; B := [31, 27, 38]; Cost := [23, 15, 6, 7, 1, 53, 4]; NonNegSol(A, B); ScalarProduct(It, Cost); Alias IP := $contrib/intprog; IP.Sol(A, B, Cost); TS := IP.TestSet(A, Cost); IP.Sol(TS, B, Cost); IP.Sol(TS, [28, 23, 33], Cost); Comp(IP.Sol(TS, B, Cost), 'MinSol'); Comp(IP.Sol(TS, B, Cost), 'MinCost'); -- You can directly use the GBasis/Test-Set in its ring: TS := IP.TestSet(A, Cost); Use Var(RingEnv(TS.TestSet[1])); TS; NR(LogToTerm(NonNegSol(A, B)), TS.TestSet); -- If A contains negative entries, expect a much slower computation -- it will not use the optimized algorithm 'Toric' A := Mat([ [ 1, 2, 3, 4], [-1,-1, 0, 1]]); B := [30,-10]; Cost := [ 1, 1, 1, 1]; IP.Sol(A, B, Cost); IP.SolIneq(A, B, Cost); " EndDefine; -- Man ------[ Main functions ]-------- Define TestSet(X, Cost) If Type(X) = MAT Then X := IP.Tagged(Record(A=X)) ElsIf Type(X) <> TAGGED(IP.Tag()) Then Error("TestSet: first argument must be MAT or ", IP.Tag()) EndIf; Return IP.AuxTestSet(X, Cost, 0) EndDefine; -- TestSet Define Sol(X, B, Cost) If Type(X) = MAT Then If Min(Flatten(List(X)))<0 Then Return IP.SolNeg(X,B,Cost) Else X := IP.Tagged(Record(A=X)) EndIf; ElsIf Type(X) <> TAGGED(IP.Tag()) Then Error("Sol: first argument must be MAT or ", IP.Tag()) EndIf; IPSol := NonNegSol(X.A, B); If IPSol=[] Then Return Record(MinSol=[], MinCost=NULL) EndIf; X := IP.AuxTestSet(X, Cost, Deg(LogToTerm(B))); MinSol := IP.NF(X, IPSol); Return Record(MinSol = MinSol, MinCost = ScalarProduct(MinSol, X.Cost)) EndDefine; -- Sol Define SolIneq(A, B, Cost) I := Identity(Len(A)); AI := Mat([ Concat(A[J], I[J]) | J In 1..Len(A) ]); SolEq := IP.SolNeg(AI, B, Concat(Cost, NewList(Len(A), 0))); MinSol := First(SolEq.MinSol, Len(A[1])); Return Record(MinSol = MinSol, MinCost = ScalarProduct(MinSol, Cost)) EndDefine; -- SolIneq ------[ pretty printing ]-------- Tag() := IP.PkgName() + ".IP"; Tagged(X) := Tagged(X, IP.Tag()); Define Print_IP(IP_Rec) PrintLn "IP.Tagged( Record["; F := Fields(IP_Rec); Foreach I In 1..(Len(F)-1) Do PrintLn " ", F[I], " = ", IP_Rec.Var(F[I]), "," EndForeach; /* TestSet */ R := RingEnv(Comp(IP_Rec.Var(F[Len(F)]), 1)); PrintLn " ", F[Len(F)], " = ", R, "::"; Using Var(R) Do PrintLn IP_Rec.Var(F[Len(F)]); EndUsing; PrintLn " ])" EndDefine; -- Print_IP ------[ Auxiliary functions ]-------- Define AuxTestSet(IP, Cost, Trunc) If (IsDefined(IP.Cost) And IP.Cost=Cost And IsDefined(IP.TestSet)) And ( Not IsDefined(IP.Trunc) Or IP.Trunc>=Trunc ) Then Return IP EndIf; IP.Cost := Cost; IPRingName := IP.SuggestedRing(IP.A, IP.Cost); If Not IsDefined(IP.Ideal) Then I := Var(IPRingName) :: Toric(IP.A); IP.Ideal := I; ElsIf RingEnv(IP.Ideal)<>IPRingName Then I := Var(IPRingName) :: Image(IP.Ideal.Gens, x); Else I := IP.Ideal; EndIf; IP.TestSet := Var(IPRingName) :: TestSet(I.Gens, Trunc); Return IP EndDefine; -- AuxTestSet Define VectToBin(V) Return LogToTerm([Max(0,V[I]) | I In 1..Len(V)]) - LogToTerm([-Min(0,V[I]) | I In 1..Len(V)]); EndDefine; -- VectToBin Define HomogWeights(A) W := Sum([ R | R In A ]); If W = NewList(Len(W), W[1]) Then Return NewList(Len(W), 1) Else Foreach R In A Do If Product(R)<>0 Then Return R EndIf; EndForeach; EndIf; Return W EndDefine; -- HomogWeights Define Ord(Ws, Cost) LenWs := Len(Ws); O := NewMat(LenWs, LenWs, 0); O[1] := NewList(LenWs, 1); O[2] := [Cost[I]/Ws[I] | I In 1.. LenWs]; DenLCM := 1; For I := 1 To LenWs Do If Type(O[2,I])=RAT Then DenLCM := LCM(DenLCM, O[2,I].Den) EndIf; EndFor; If DenLCM<>1 Then For I := 1 To LenWs Do O[2,I] := DenLCM*O[2,I] EndFor; EndIf; For I := 3 To LenWs Do O[I, LenWs+3-I] := -1 EndFor; If O[2,1]<>O[2,2] Then Return O EndIf; I := 3; Row := 2; While I<=LenWs And Row=2 Do If O[2,1]<>O[2,I] Then Row:=LenWs+3-I EndIf; I := I+1; EndWhile; For I := Row To LenWs-1 Do O[I] := O[I+1] EndFor; O[LenWs, 2] := -1; O[LenWs, 3] := 0; Return O EndDefine; -- Ord Define SuggestedRing(A, Cost) Ws := IP.HomogWeights(A); O := IP.Ord(Ws, Cost); IPRingName := "IPRing" + Sum(["_"+Sprint(X)|X In Ws]) + Sum(["_"+Sprint(X)|X In Cost]); If Not(IPRingName IsIn RingEnvs() ) Then Var(IPRingName) ::= CoeffRing[x[1..Len(Ws)]], Weights(Ws), Ord(O) EndIf; Return IPRingName EndDefine; -- SuggestedRing -------------------------------------------------------------- Define OrdMat2(A, Cost) R := Len(A); C := Len(A[1]); M := BlockMatrix([ [[NewList(R+1,1)], 0], [ 0, [Cost]], [First(LexMat(R+1), R), 0], [0, Last(LexMat(C), C-1)]]); If Det(M) <> 0 Then Return M EndIf; I := Min([J In 1..Len(Cost) | Cost[J]<>0]); MM := BlockMatrix([ [First(LexMat(C), I-1)], [Last(LexMat(C), C-I)]]); Return BlockMatrix([ [[NewList(R+1,1)], 0], [ 0, [Cost]], [First(LexMat(R+1), R), 0], [0, MM]]); EndDefine; -- OrdMat2 Define SolNeg(A, B, Cost) If Min(Cost)<0 Then Error("Sol: matrix has negative entries; cost function must be positive") EndIf; MinB := Min(B); If MinB > 0 Then MinB := 0 EndIf; Ring := NewId(); Var(Ring) ::= Z/(2)[ z[1..Len(A)] t y[1..Len(A[1])] ], Ord(IP.OrdMat2(A,Cost)); Using Var(Ring) Do VB := t^(-MinB) * Product([ z[J]^(B[J]-MinB) | J In 1..Len(A)]); F := NF(VB, IP.Ideal2(A)); If LT(F) < t Then MinSol := Last(Log(F), Len(A[1])); Return Record(MinSol = MinSol, MinCost = ScalarProduct(MinSol, Cost)) EndIf; EndUsing; Return []; EndDefine; -- SolNeg Define Ideal2(A) C := Len(A[1]); R := Len(A); L := [ Product([z[I]^Max( A[I,J],0) | I In 1..R]) - y[J]*Product([z[I]^Max(-A[I,J],0) | I In 1..R]) | J In 1..C ]; Append(L, 1-t * Product([z[I] | I In 1..R])); Return Ideal(L); EndDefine; -- Ideal2 -------------------------------------------------------------- Define NF(X,IPSol) Return Var(RingEnv(X.TestSet[1]))::Log(NR(LogToTerm(IPSol),X.TestSet)); EndDefine; -- NF ------[ CheckInput ]-------- EndPackage; -- Package $contrib/intprog