login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit

)abbrev package DOPT DiscreteOptimalTransport
++ Author: Kurt Pagani
++ Date Created: Thu Jun 21 01:42:00 CEST 2018
++ References:
++ Description:
++
DiscreteOptimalTransport(R) : Exports == Implementation where
R:Join(Ring,IntegralDomain,OrderedSet)
PI ==> PositiveInteger
NNI ==> NonNegativeInteger
INT ==> Integer
LLI ==> List List Integer
MATR ==> Matrix R
VECR ==> Vector R
null ==> empty?
Exports ==  with
constraintMatrix : (PI,PI) -> MATR
northWest : (PI,PI,VECR,VECR) -> MATR
++ NW rule
getBasis : MATR -> List List List Integer
projRows : (LLI,INT) -> List Integer
++ Given a basis B and an integer c, return a list of integers
++ [r1,..,rn], such that each (ri,c) is in B, i=1..n.
projCols : (LLI,INT) -> List Integer
++ Given a basis B and an integer r, return a list of integers
++ [c1,..,cn], such that each (r,ci) is in B, i=1..n.
recSolve : (PI,PI,LLI,MATR) -> List VECR
computeReducedCosts : (PI,PI,List VECR,LLI,MATR)-> MATR
determineAlpha : (LLI,MATR,MATR) -> Record(alpha:R,ijmax:List INT)
checkHoriz : (LLI,MATR,LLI,INT,INT,INT,INT) -> LLI
checkVert  : (LLI,MATR,LLI,INT,INT,INT,INT) -> LLI
findPath : (MATR,LLI,INT,INT) -> LLI
findDelta : (LLI,MATR) -> Record(delta0:R,ij0:List INT)
computeSolutionMatrix : (PI,PI,VECR,VECR,MATR) -> MATR
modiMethod :(VECR,VECR,MATR) -> Record(Sol:MATR,cost:R)
--coerce : % -> OutputForm
constraintMatrix(m:PI,n:PI):MATR ==
M:MATR:=new(m,n,0$R)$MATR
e:VECR:=new(n,1$R)$VECR
LM:List MATR:=[setRow!(copy M,i,e) for i in 1..m]
U:MATR:=horizConcat(LM)$MATR E:MATR:=diagonalMatrix(e) LD:List MATR:=[E for i in 1..m] D:=horizConcat(LD)$MATR
vertConcat(U,D)
northWest(m:PI,n:PI,a:VECR,b:VECR):MATR ==
X:MATR:=new(m,n,0$R)$MATR
aa:VECR:=new(n,0$R)$VECR
bb:VECR:=new(m,0$R)$VECR
u:INT:=1 ; v:INT:=1
while u<=m and v<=n repeat
if b.v - aa.v < a.u - bb.u then
z:=b.v-aa.v
X(u,v):=z
aa.v := aa.v + z
bb.u := bb.u + z
v := v+ 1
else
z:=a.u-bb.u
X(u,v):=z
aa.v := aa.v + z
bb.u := bb.u + z
u := u+ 1
return X
getBasis(X:MATR):List List List Integer ==
b:List List Integer:=[]
c:List List Integer:=[]
for i in 1..nrows(X) repeat
for j in 1..ncols(X) repeat
if X(i,j) = 0$R then c:=concat(c,[[i,j]]) else b:=concat(b,[[i,j]]) return [b,c] projRows(b:LLI,nc:INT):List Integer == [first(t) for t in b | t.2=nc] projCols(b:LLI,nr:INT):List Integer == [second(t) for t in b | t.1=nr] recSolve(m:PI,n:PI,b:LLI,C:MATR):List VECR == u:VECR:=new(m,0$R)$VECR v:VECR:=new(n,0$R)$VECR ud:List INT:=[i for i in 1..m] vd:List INT:=[i for i in 1..n] bn:List Integer:=projRows(b,n) for i in bn repeat u.i:=C(i,n) ud:=remove!(i,ud) for j in projCols(b,i) repeat v.j:=C(i,j)-u.i vd:=remove!(j,vd) repeat for k in ud repeat for j in projCols(b,k) repeat if not member?(j,vd) then u.k:=C(k,j)-v.j ud:=remove!(k,ud) for l in vd repeat for j in projRows(b,l) repeat if not member?(j,ud) then v.l:=C(j,l)-u.j vd:=remove!(l,vd) if null(ud) and null(vd) then return [u,v] computeReducedCosts(m:PI,n:PI,rs:List VECR,c:LLI,C:MATR):MATR == u:VECR:=rs.1 v:VECR:=rs.2 RC:MATR:=new(m,n,0$R)$MATR for t in c repeat RC(t.1,t.2):=u.(t.1)+v.(t.2) return RC determineAlpha(c:LLI,RC:MATR,C:MATR):Record(alpha:R,ijmax:List INT) == mm:R:=RC(c.1.1,c.1.2)-C(c.1.1,c.1.2) tmax:List Integer:=c.1 for t in c repeat x:R:=RC(t.1,t.2)-C(t.1,t.2) if x>mm then mm:=x tmax:=t return [mm,tmax] checkHoriz(path:LLI,X:MATR,bb:LLI,u:INT,v:INT,u1:INT,v1:INT):LLI == n:=ncols X for i in 1..n repeat if i ~= v and member?([u,i],bb) then if i=v1 then path:=concat(path,[u,i]) return path p := checkVert(path,X,bb,u,i,u1,v1) if not null p then path:=concat(p,[u,i]) return path return []::LLI checkVert(path:LLI,X:MATR,bb:LLI,u:INT,v:INT,u1:INT,v1:INT):LLI == m:=nrows X for i in 1..m repeat if i ~= u and member?([i,v],bb) then p := checkHoriz(path,X,bb,i,v,u1,v1) if not null p then path:=concat(p,[i,v]) return path return []::LLI findPath(X:MATR,bb:LLI,u:INT,v:INT):LLI == p0:LLI:=[[u,v]] path:LLI:=checkHoriz(p0,X,bb,u,v,u,v) return path findDelta(cy:LLI,X:MATR):Record(delta0:R,ij0:List INT) == ncy:LLI:=[cy.j for j in 1..#cy | even? j] --nX:List R:=[X(t.1,t.2) for t in ncy] tmin:List Integer:=ncy.1 m:R:=X(tmin.1,tmin.2) for t in ncy repeat x:R:=X(t.1,t.2) if x<m then m:=x tmin:=t return [m,tmin] maxIter:INT:=10 computeSolutionMatrix(m:PI,n:PI,a:VECR,b:VECR,C:MATR):MATR == X:MATR:=northWest(m,n,a,b) bl:List LLI:=getBasis(X) bb:LLI:=bl.1 cc:LLI:=bl.2 it:INT:=0 -- repeat it:=it+1 if it > maxIter then break rs:List VECR:=recSolve(m,n,bb,C) RC:MATR:=computeReducedCosts(m,n,rs,cc,C) -- da:=determineAlpha(cc,RC,C) -- da.alpha, da.ijmax if da.alpha > 0$R then
path:=findPath(X,bb,da.ijmax.1,da.ijmax.2)
d0:=findDelta(path,X) -- d0.delta0, d0.ij0
for j in 1..#path repeat
p:=path.j
if odd?(j) then
X(p.1,p.2):=X(p.1,p.2)+d0.delta0
else
X(p.1,p.2):=X(p.1,p.2)-d0.delta0
bb:=remove(d0.ij0,bb)
bb:=concat(bb,da.ijmax)
cc:=remove(da.ijmax,cc)
cc:=concat(cc,d0.ij0)
else
return X
X
modiMethod(a:VECR,b:VECR,C:MATR):Record(Sol:MATR,cost:R) ==
m0:NNI:=nrows(C)
n0:NNI:=ncols(C)
if m0>0 then m:PI:=m0::PI else error "1"
if n0>0 then n:PI:=n0::PI else error "2"
if #a~=m then error "3"
if #b~=n then error "4"
if reduce(_+,a) ~= reduce(_+,b) then error "5"
S:MATR:=computeSolutionMatrix(m,n,a,b,C)
CTS:MATR:=transpose(C)*S
c:R:=reduce(_+,[CTS(i,i) for i in 1..n])  -- trace reqs too much
return [S,c]
   Compiling FriCAS source code from file
using old system compiler.
DOPT abbreviates package DiscreteOptimalTransport
------------------------------------------------------------------------
initializing NRLIB DOPT for DiscreteOptimalTransport
compiling into NRLIB DOPT
****** Domain: R already in scope
****** Domain: R already in scope
compiling exported constraintMatrix : (PositiveInteger,PositiveInteger) -> Matrix R
Time: 0.04 SEC.
compiling exported northWest : (PositiveInteger,PositiveInteger,Vector R,Vector R) -> Matrix R
Time: 0.02 SEC.
compiling exported getBasis : Matrix R -> List List List Integer
Time: 0.01 SEC.
compiling exported projRows : (List List Integer,Integer) -> List Integer
Time: 0.01 SEC.
compiling exported projCols : (List List Integer,Integer) -> List Integer
Time: 0 SEC.
compiling exported recSolve : (PositiveInteger,PositiveInteger,List List Integer,Matrix R) -> List Vector R
Time: 0.04 SEC.
compiling exported computeReducedCosts : (PositiveInteger,PositiveInteger,List Vector R,List List Integer,Matrix R) -> Matrix R
Time: 0.04 SEC.
compiling exported determineAlpha : (List List Integer,Matrix R,Matrix R) -> Record(alpha: R,ijmax: List Integer)
Time: 0.01 SEC.
compiling exported checkHoriz : (List List Integer,Matrix R,List List Integer,Integer,Integer,Integer,Integer) -> List List Integer
Time: 0.01 SEC.
compiling exported checkVert : (List List Integer,Matrix R,List List Integer,Integer,Integer,Integer,Integer) -> List List Integer
Time: 0.01 SEC.
compiling exported findPath : (Matrix R,List List Integer,Integer,Integer) -> List List Integer
Time: 0.01 SEC.
compiling exported findDelta : (List List Integer,Matrix R) -> Record(delta0: R,ij0: List Integer)
Time: 0.01 SEC.
compiling exported computeSolutionMatrix : (PositiveInteger,PositiveInteger,Vector R,Vector R,Matrix R) -> Matrix R
Time: 0.11 SEC.
compiling exported modiMethod : (Vector R,Vector R,Matrix R) -> Record(Sol: Matrix R,cost: R)
Time: 0.02 SEC.
(time taken in buildFunctor:  0)
;;;     ***       |DiscreteOptimalTransport| REDEFINED
;;;     ***       |DiscreteOptimalTransport| REDEFINED
Time: 0 SEC.
Warnings:
[1] getBasis:  b has no value
[2] getBasis:  c has no value
[3] recSolve:  ud has no value
[4] recSolve:  vd has no value
[5] determineAlpha:  mm has no value
[6] determineAlpha:  tmax has no value
[7] findDelta:  m has no value
[8] findDelta:  tmin has no value
[9] computeSolutionMatrix:  alpha has no value
[10] computeSolutionMatrix:  ijmax has no value
[11] computeSolutionMatrix:  delta0 has no value
[12] computeSolutionMatrix:  ij0 has no value
Cumulative Statistics for Constructor DiscreteOptimalTransport
Time: 0.34 seconds
finalizing NRLIB DOPT
Processing DiscreteOptimalTransport for Browser database:
--------constructor---------
--->-->DiscreteOptimalTransport((constraintMatrix ((Matrix R) (PositiveInteger) (PositiveInteger)))): Not documented!!!!
--------(northWest ((Matrix R) (PositiveInteger) (PositiveInteger) (Vector R) (Vector R)))---------
--->-->DiscreteOptimalTransport((getBasis ((List (List (List (Integer)))) (Matrix R)))): Not documented!!!!
--------(projRows ((List (Integer)) (List (List (Integer))) (Integer)))---------
--->-->DiscreteOptimalTransport((projRows ((List (Integer)) (List (List (Integer))) (Integer)))): Improper first word in comments: Given
--------(projCols ((List (Integer)) (List (List (Integer))) (Integer)))---------
--->-->DiscreteOptimalTransport((projCols ((List (Integer)) (List (List (Integer))) (Integer)))): Improper first word in comments: Given
--->-->DiscreteOptimalTransport((recSolve ((List (Vector R)) (PositiveInteger) (PositiveInteger) (List (List (Integer))) (Matrix R)))): Not documented!!!!
--->-->DiscreteOptimalTransport((computeReducedCosts ((Matrix R) (PositiveInteger) (PositiveInteger) (List (Vector R)) (List (List (Integer))) (Matrix R)))): Not documented!!!!
--->-->DiscreteOptimalTransport((determineAlpha ((Record (: alpha R) (: ijmax (List (Integer)))) (List (List (Integer))) (Matrix R) (Matrix R)))): Not documented!!!!
--->-->DiscreteOptimalTransport((checkHoriz ((List (List (Integer))) (List (List (Integer))) (Matrix R) (List (List (Integer))) (Integer) (Integer) (Integer) (Integer)))): Not documented!!!!
--->-->DiscreteOptimalTransport((checkVert ((List (List (Integer))) (List (List (Integer))) (Matrix R) (List (List (Integer))) (Integer) (Integer) (Integer) (Integer)))): Not documented!!!!
--->-->DiscreteOptimalTransport((findPath ((List (List (Integer))) (Matrix R) (List (List (Integer))) (Integer) (Integer)))): Not documented!!!!
--->-->DiscreteOptimalTransport((findDelta ((Record (: delta0 R) (: ij0 (List (Integer)))) (List (List (Integer))) (Matrix R)))): Not documented!!!!
--->-->DiscreteOptimalTransport((computeSolutionMatrix ((Matrix R) (PositiveInteger) (PositiveInteger) (Vector R) (Vector R) (Matrix R)))): Not documented!!!!
--->-->DiscreteOptimalTransport((modiMethod ((Record (: Sol (Matrix R)) (: cost R)) (Vector R) (Vector R) (Matrix R)))): Not documented!!!!
; compiling file "/var/aw/var/LatexWiki/DOPT.NRLIB/DOPT.lsp" (written 16 JUL 2018 07:46:12 PM):
; /var/aw/var/LatexWiki/DOPT.NRLIB/DOPT.fasl written
; compilation finished in 0:00:00.177
------------------------------------------------------------------------
DiscreteOptimalTransport is now explicitly exposed in frame initial
DiscreteOptimalTransport will be automatically loaded when needed
from /var/aw/var/LatexWiki/DOPT.NRLIB/DOPT

Test flavours

fricas
--)co dopt
supply, demand a:=vector [4,7,4]$Vector INT  (1) Type: Vector(Integer) fricas b:=vector [2,4,6,3]$Vector INT
 (2)
Type: Vector(Integer)
fricas
-- cost matrix
C:= matrix [[2,3,4,1],[5,4,2,3],[4,2,8,6]]
 (3)
Type: Matrix(Integer)
fricas
X:=computeSolutionMatrix(3,4,a,b,C)
 (4)
Type: Matrix(Integer)
fricas
cost:=trace(transpose(C)*X) -- 29
 (5)
Type: PositiveInteger?
fricas
--2
C2 := matrix [[ 2, 1, 3, 3, 2, 5],[ 3, 2, 2, 4, 3, 4], _
[ 3, 5, 4, 2, 4, 1],[ 4, 2, 2, 1, 2, 2]]
 (6)
Type: Matrix(Integer)
fricas
b2:= vector [ 30, 50, 20, 40, 30, 11]
 (7)
Type: Vector(PositiveInteger?)
fricas
a2:= vector [ 50, 40, 60, 31]
 (8)
Type: Vector(PositiveInteger?)
fricas
X2:=computeSolutionMatrix(4,6,a2,b2,C2)
 (9)
Type: Matrix(Integer)
fricas
cost2:=trace(transpose(C2)*X2) -- 330
 (10)
Type: PositiveInteger?
fricas
-- 3
C3:= matrix [[ 1, 2, 1, 4, 5, 2],[ 3, 3, 2, 1, 4, 3], _
[ 4, 2, 5, 9, 6, 2],[ 3, 1, 7, 3, 4, 6]]
 (11)
Type: Matrix(Integer)
fricas
b3:= vector [ 20, 40, 30, 10, 50, 25]
 (12)
Type: Vector(PositiveInteger?)
fricas
a3:= vector [ 30, 50, 75, 20]
 (13)
Type: Vector(PositiveInteger?)
fricas
X3:=computeSolutionMatrix(4,6,a3,b3,C3)
 (14)
Type: Matrix(Integer)
fricas
cost3:=trace(transpose(C3)*X3) -- 430
 (15)
Type: PositiveInteger?
fricas
-- 4 ??? sum(a)~sum(b)
C4:= matrix [[ 5, 3, 6, 2],[ 4, 7, 9, 1],[ 3, 4, 7, 5]]
 (16)
Type: Matrix(Integer)
fricas
b4 := vector [ 16, 18, 30, 25]
 (17)
Type: Vector(PositiveInteger?)
fricas
a4 := vector [ 19, 37, 34]
 (18)
Type: Vector(PositiveInteger?)
fricas
X4:=computeSolutionMatrix(3,4,a4,b4,C4)
 (19)
Type: Matrix(Integer)
fricas
cost4:=trace(transpose(C4)*X4) -- 348
 (20)
Type: PositiveInteger?
fricas
-- 5 WP TT
C5:= matrix [[6,4,3],[3,2,2]]
 (21)
Type: Matrix(Integer)
fricas
a5:= vector [16,14]
 (22)
Type: Vector(PositiveInteger?)
fricas
b5:= vector [15,5,10]
 (23)
Type: Vector(PositiveInteger?)
fricas
X5:=computeSolutionMatrix(2,3,a5,b5,C5)
 (24)
Type: Matrix(Integer)
fricas
cost5:=trace(transpose(C5)*X5) -- 98
 (25)
Type: PositiveInteger?
fricas
-- 6 WP ZykM
C6:= matrix [[7,4,3],[5,5,6]]
 (26)
Type: Matrix(Integer)
fricas
a6:= vector [12,8]
 (27)
Type: Vector(PositiveInteger?)
fricas
b6:= vector [4,10,6]
 (28)
Type: Vector(PositiveInteger?)
fricas
X6:=computeSolutionMatrix(2,3,a6,b6,C6)
 (29)
Type: Matrix(Integer)
fricas
cost6:=trace(transpose(C6)*X6) -- 82
 (30)
Type: PositiveInteger?
fricas
--- 7
C7:= matrix [[140,130,170,150,140],[150,160,140,180,130],[160,120,180,170,190]]
 (31)
Type: Matrix(Integer)
fricas
a7:= vector [900,800,400]
 (32)
Type: Vector(PositiveInteger?)
fricas
b7:= vector [300,500,400,600,300]
 (33)
Type: Vector(PositiveInteger?)
fricas
X7:=computeSolutionMatrix(3,5,a7,b7,C7)
 (34)
Type: Matrix(Integer)
fricas
cost7:=trace(transpose(C7)*X7) -- 289000
 (35)
Type: PositiveInteger?
fricas
--- 8
C8:= matrix [[3,2],[1,5],[5,4]]
 (36)
Type: Matrix(Integer)
fricas
a8:= vector [45,60,5]
 (37)
Type: Vector(PositiveInteger?)
fricas
b8:= vector [50,60]
 (38)
Type: Vector(PositiveInteger?)
fricas
X8:=computeSolutionMatrix(3,2,a8,b8,C8)
 (39)
Type: Matrix(Integer)
fricas
cost8:=trace(transpose(C8)*X8)  -- 210
 (40)
Type: PositiveInteger?
fricas
--- 9
C9:= matrix [[5,5,2,7],[8,3,1,8],[2,4,4,5]]
 (41)
Type: Matrix(Integer)
fricas
a9:= vector [9,15,14]
 (42)
Type: Vector(PositiveInteger?)
fricas
b9:= vector [10,12,9,7]
 (43)
Type: Vector(PositiveInteger?)
fricas
X9:=computeSolutionMatrix(3,4,a9,b9,C9)
 (44)
Type: Matrix(Integer)
fricas
cost9:=trace(transpose(C9)*X9)  -- 112
 (45)
Type: PositiveInteger?
fricas
modiMethod(a9,b9,C9) -- check sum a = sum b
 (46)
Type: Record(Sol: Matrix(Integer),cost: Integer)

Note: "null" must have been removed recently ?? WA: null ==> empty?

 Subject:   Be Bold !! ( 14 subscribers )