help links subscribe changes refresh edit

	SandBox
SandBoxDGE ←	↑	→ SandBoxDefineInteger

SandBoxDOPT

Edit detail for SandBoxDOPT revision 1 of 1

	1
Editor: pagani Time: 2018/07/16 19:46:12 GMT+0
Note:

changed:
-
\begin{spad}
)abbrev package DOPT DiscreteOptimalTransport
++ Author: Kurt Pagani
++ Date Created: Thu Jun 21 01:42:00 CEST 2018
++ License: BSD
++ 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
	
	
  Implementation ==  add 
    
    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]
      

\end{spad}

Test flavours

\begin{axiom}
--)co dopt

-- supply, demand
a:=vector [4,7,4]$Vector INT
b:=vector [2,4,6,3]$Vector INT
-- cost matrix
C:= matrix [[2,3,4,1],[5,4,2,3],[4,2,8,6]]

X:=computeSolutionMatrix(3,4,a,b,C)
cost:=trace(transpose(C)*X) -- 29

--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]]

b2:= vector [ 30, 50, 20, 40, 30, 11]
a2:= vector [ 50, 40, 60, 31]

X2:=computeSolutionMatrix(4,6,a2,b2,C2)
cost2:=trace(transpose(C2)*X2) -- 330

-- 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]]
             
b3:= vector [ 20, 40, 30, 10, 50, 25]
a3:= vector [ 30, 50, 75, 20]

X3:=computeSolutionMatrix(4,6,a3,b3,C3)
cost3:=trace(transpose(C3)*X3) -- 430

-- 4 ??? sum(a)~sum(b)

C4:= matrix [[ 5, 3, 6, 2],[ 4, 7, 9, 1],[ 3, 4, 7, 5]]
b4 := vector [ 16, 18, 30, 25]
a4 := vector [ 19, 37, 34] 

X4:=computeSolutionMatrix(3,4,a4,b4,C4)
cost4:=trace(transpose(C4)*X4) -- 348

-- 5 WP TT
C5:= matrix [[6,4,3],[3,2,2]]
a5:= vector [16,14]
b5:= vector [15,5,10]

X5:=computeSolutionMatrix(2,3,a5,b5,C5)
cost5:=trace(transpose(C5)*X5) -- 98

-- 6 WP ZykM
C6:= matrix [[7,4,3],[5,5,6]]
a6:= vector [12,8]
b6:= vector [4,10,6]

X6:=computeSolutionMatrix(2,3,a6,b6,C6)
cost6:=trace(transpose(C6)*X6) -- 82

--- 7

C7:= matrix [[140,130,170,150,140],[150,160,140,180,130],[160,120,180,170,190]]
a7:= vector [900,800,400]
b7:= vector [300,500,400,600,300]

X7:=computeSolutionMatrix(3,5,a7,b7,C7)
cost7:=trace(transpose(C7)*X7) -- 289000

--- 8
C8:= matrix [[3,2],[1,5],[5,4]]
a8:= vector [45,60,5]
b8:= vector [50,60]

X8:=computeSolutionMatrix(3,2,a8,b8,C8)
cost8:=trace(transpose(C8)*X8)  -- 210

--- 9
C9:= matrix [[5,5,2,7],[8,3,1,8],[2,4,4,5]]
a9:= vector [9,15,14]
b9:= vector [10,12,9,7]

X9:=computeSolutionMatrix(3,4,a9,b9,C9)
cost9:=trace(transpose(C9)*X9)  -- 112
modiMethod(a9,b9,C9) -- check sum a = sum b


-- modiMethod(a,b,C) ... modiMethod(a8,b8,C8) 
\end{axiom}

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

fricas

(1) -> <spad>

fricas

)abbrev package DOPT DiscreteOptimalTransport
++ Author: Kurt Pagani
++ Date Created: Thu Jun 21 01:42:00 CEST 2018
++ License: BSD
++ 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

  Implementation ==  add 

    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]</spad>

fricas

Compiling FriCAS source code from file 
      /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/4555332396781246409-25px001.spad
      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.03 SEC.

   compiling exported northWest : (PositiveInteger,PositiveInteger,Vector R,Vector R) -> Matrix R
Time: 0.01 SEC.

   compiling exported getBasis : Matrix R -> List List List Integer
Time: 0 SEC.

   compiling exported projRows : (List List Integer,Integer) -> List Integer
Time: 0 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.02 SEC.

   compiling exported computeReducedCosts : (PositiveInteger,PositiveInteger,List Vector R,List List Integer,Matrix R) -> Matrix R
Time: 0.02 SEC.

   compiling exported determineAlpha : (List List Integer,Matrix R,Matrix R) -> Record(alpha: R,ijmax: List Integer)
Time: 0 SEC.

   compiling exported checkHoriz : (List List Integer,Matrix R,List List Integer,Integer,Integer,Integer,Integer) -> List List Integer
Time: 0 SEC.

   compiling exported checkVert : (List List Integer,Matrix R,List List Integer,Integer,Integer,Integer,Integer) -> List List Integer
Time: 0 SEC.

   compiling exported findPath : (Matrix R,List List Integer,Integer,Integer) -> List List Integer
Time: 0 SEC.

   compiling exported findDelta : (List List Integer,Matrix R) -> Record(delta0: R,ij0: List Integer)
Time: 0 SEC.

   compiling exported computeSolutionMatrix : (PositiveInteger,PositiveInteger,Vector R,Vector R,Matrix R) -> Matrix R
Time: 0.04 SEC.

   compiling exported modiMethod : (Vector R,Vector R,Matrix R) -> Record(Sol: Matrix R,cost: R)
Time: 0 SEC.

(time taken in buildFunctor:  0)
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

   Cumulative Statistics for Constructor DiscreteOptimalTransport
      Time: 0.17 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
"Given a basis \\spad{B} and an integer \\spad{c},{} return a list of integers [\\spad{r1},{}..,{}\\spad{rn}],{} such that each (\\spad{ri},{}\\spad{c}) is in \\spad{B},{} \\spad{i=1}..\\spad{n}."
--------(projCols ((List (Integer)) (List (List (Integer))) (Integer)))---------
--->-->DiscreteOptimalTransport((projCols ((List (Integer)) (List (List (Integer))) (Integer)))): Improper first word in comments: Given
"Given a basis \\spad{B} and an integer \\spad{r},{} return a list of integers [\\spad{c1},{}..,{}\\spad{cn}],{} such that each (\\spad{r},{}\\spad{ci}) is in \\spad{B},{} \\spad{i=1}..\\spad{n}."
--->-->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 13 OCT 2025 07:21:01 AM):

; wrote /var/aw/var/LatexWiki/DOPT.NRLIB/DOPT.fasl
; compilation finished in 0:00:00.132
------------------------------------------------------------------------
   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

$\label{eq1}\left[ 4, \: 7, \: 4 \right]$

(1)

Type: Vector(Integer)

fricas

b:=vector [2,4,6,3]$Vector INT

$\label{eq2}\left[ 2, \: 4, \: 6, \: 3 \right]$

(2)

Type: Vector(Integer)

fricas

-- cost matrix
C:= matrix [[2,3,4,1],[5,4,2,3],[4,2,8,6]]

$\label{eq3}\left[ \begin{array}{cccc} 2 & 3 & 4 & 1 \ 5 & 4 & 2 & 3 \ 4 & 2 & 8 & 6$

(3)

Type: Matrix(Integer)

fricas

X:=computeSolutionMatrix(3,4,a,b,C)

$\label{eq4}\left[ \begin{array}{cccc} 2 & 0 & 0 & 2 \ 0 & 0 & 6 & 1 \ 0 & 4 & 0 & 0$

(4)

Type: Matrix(Integer)

fricas

cost:=trace(transpose(C)*X) -- 29

$\label{eq5}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]]

$\label{eq6}\left[ \begin{array}{cccccc} 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]

$\label{eq7}\left[{30}, \:{50}, \:{20}, \:{40}, \:{30}, \:{11}\right]$

(7)

Type: Vector(PositiveInteger?)

fricas

a2:= vector [ 50, 40, 60, 31]

$\label{eq8}\left[{50}, \:{40}, \:{60}, \:{31}\right]$

(8)

Type: Vector(PositiveInteger?)

fricas

X2:=computeSolutionMatrix(4,6,a2,b2,C2)

$\label{eq9}\left[ \begin{array}{cccccc} {20}&{30}& 0 & 0 & 0 & 0 \ 0 &{20}&{20}& 0 & 0 & 0 \ {10}& 0 & 0 &{39}& 0 &{11} \ 0 & 0 & 0 & 1 &{30}& 0$

(9)

Type: Matrix(Integer)

fricas

cost2:=trace(transpose(C2)*X2) -- 330

$\label{eq10}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]]

$\label{eq11}\left[ \begin{array}{cccccc} 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]

$\label{eq12}\left[{20}, \:{40}, \:{30}, \:{10}, \:{50}, \:{25}\right]$

(12)

Type: Vector(PositiveInteger?)

fricas

a3:= vector [ 30, 50, 75, 20]

$\label{eq13}\left[{30}, \:{50}, \:{75}, \:{20}\right]$

(13)

Type: Vector(PositiveInteger?)

fricas

X3:=computeSolutionMatrix(4,6,a3,b3,C3)

$\label{eq14}\left[ \begin{array}{cccccc} {20}& 0 &{10}& 0 & 0 & 0 \ 0 & 0 &{20}&{10}&{20}& 0 \ 0 &{40}& 0 & 0 &{10}&{25} \ 0 & 0 & 0 & 0 &{20}& 0$

(14)

Type: Matrix(Integer)

fricas

cost3:=trace(transpose(C3)*X3) -- 430

$\label{eq15}430$

(15)

Type: PositiveInteger?

fricas

-- 4 ??? sum(a)~sum(b)

C4:= matrix [[ 5, 3, 6, 2],[ 4, 7, 9, 1],[ 3, 4, 7, 5]]

$\label{eq16}\left[ \begin{array}{cccc} 5 & 3 & 6 & 2 \ 4 & 7 & 9 & 1 \ 3 & 4 & 7 & 5$

(16)

Type: Matrix(Integer)

fricas

b4 := vector [ 16, 18, 30, 25]

$\label{eq17}\left[{16}, \:{18}, \:{30}, \:{25}\right]$

(17)

Type: Vector(PositiveInteger?)

fricas

a4 := vector [ 19, 37, 34]

$\label{eq18}\left[{19}, \:{37}, \:{34}\right]$

(18)

Type: Vector(PositiveInteger?)

fricas

X4:=computeSolutionMatrix(3,4,a4,b4,C4)

$\label{eq19}\left[ \begin{array}{cccc} 0 &{18}& 1 & 0 \ {12}& 0 & 0 &{25} \ 4 & 0 &{29}& 0$

(19)

Type: Matrix(Integer)

fricas

cost4:=trace(transpose(C4)*X4) -- 348

$\label{eq20}348$

(20)

Type: PositiveInteger?

fricas

-- 5 WP TT
C5:= matrix [[6,4,3],[3,2,2]]

$\label{eq21}\left[ \begin{array}{ccc} 6 & 4 & 3 \ 3 & 2 & 2$

(21)

Type: Matrix(Integer)

fricas

a5:= vector [16,14]

$\label{eq22}\left[{16}, \:{14}\right]$

(22)

Type: Vector(PositiveInteger?)

fricas

b5:= vector [15,5,10]

$\label{eq23}\left[{15}, \: 5, \:{10}\right]$

(23)

Type: Vector(PositiveInteger?)

fricas

X5:=computeSolutionMatrix(2,3,a5,b5,C5)

$\label{eq24}\left[ \begin{array}{ccc} 1 & 5 &{10} \ {14}& 0 & 0$

(24)

Type: Matrix(Integer)

fricas

cost5:=trace(transpose(C5)*X5) -- 98

$\label{eq25}98$

(25)

Type: PositiveInteger?

fricas

-- 6 WP ZykM
C6:= matrix [[7,4,3],[5,5,6]]

$\label{eq26}\left[ \begin{array}{ccc} 7 & 4 & 3 \ 5 & 5 & 6$

(26)

Type: Matrix(Integer)

fricas

a6:= vector [12,8]

$\label{eq27}\left[{12}, \: 8 \right]$

(27)

Type: Vector(PositiveInteger?)

fricas

b6:= vector [4,10,6]

$\label{eq28}\left[ 4, \:{10}, \: 6 \right]$

(28)

Type: Vector(PositiveInteger?)

fricas

X6:=computeSolutionMatrix(2,3,a6,b6,C6)

$\label{eq29}\left[ \begin{array}{ccc} 0 & 6 & 6 \ 4 & 4 & 0$

(29)

Type: Matrix(Integer)

fricas

cost6:=trace(transpose(C6)*X6) -- 82

$\label{eq30}82$

(30)

Type: PositiveInteger?

fricas

--- 7

C7:= matrix [[140,130,170,150,140],[150,160,140,180,130],[160,120,180,170,190]]

$\label{eq31}\left[ \begin{array}{ccccc} {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]

$\label{eq32}\left[{900}, \:{800}, \:{400}\right]$

(32)

Type: Vector(PositiveInteger?)

fricas

b7:= vector [300,500,400,600,300]

$\label{eq33}\left[{300}, \:{500}, \:{400}, \:{600}, \:{300}\right]$

(33)

Type: Vector(PositiveInteger?)

fricas

X7:=computeSolutionMatrix(3,5,a7,b7,C7)

$\label{eq34}\left[ \begin{array}{ccccc} {200}&{100}& 0 &{600}& 0 \ {100}& 0 &{400}& 0 &{300} \ 0 &{400}& 0 & 0 & 0$

(34)

Type: Matrix(Integer)

fricas

cost7:=trace(transpose(C7)*X7) -- 289000

$\label{eq35}289000$

(35)

Type: PositiveInteger?

fricas

--- 8
C8:= matrix [[3,2],[1,5],[5,4]]

$\label{eq36}\left[ \begin{array}{cc} 3 & 2 \ 1 & 5 \ 5 & 4$

(36)

Type: Matrix(Integer)

fricas

a8:= vector [45,60,5]

$\label{eq37}\left[{45}, \:{60}, \: 5 \right]$

(37)

Type: Vector(PositiveInteger?)

fricas

b8:= vector [50,60]

$\label{eq38}\left[{50}, \:{60}\right]$

(38)

Type: Vector(PositiveInteger?)

fricas

X8:=computeSolutionMatrix(3,2,a8,b8,C8)

$\label{eq39}\left[ \begin{array}{cc} 0 &{45} \ {50}&{10} \ 0 & 5$

(39)

Type: Matrix(Integer)

fricas

cost8:=trace(transpose(C8)*X8)  -- 210

$\label{eq40}210$

(40)

Type: PositiveInteger?

fricas

--- 9
C9:= matrix [[5,5,2,7],[8,3,1,8],[2,4,4,5]]

$\label{eq41}\left[ \begin{array}{cccc} 5 & 5 & 2 & 7 \ 8 & 3 & 1 & 8 \ 2 & 4 & 4 & 5$

(41)

Type: Matrix(Integer)

fricas

a9:= vector [9,15,14]

$\label{eq42}\left[ 9, \:{15}, \:{14}\right]$

(42)

Type: Vector(PositiveInteger?)

fricas

b9:= vector [10,12,9,7]

$\label{eq43}\left[{10}, \:{12}, \: 9, \: 7 \right]$

(43)

Type: Vector(PositiveInteger?)

fricas

X9:=computeSolutionMatrix(3,4,a9,b9,C9)

$\label{eq44}\left[ \begin{array}{cccc} 0 & 0 & 6 & 3 \ 0 &{12}& 3 & 0 \ {10}& 0 & 0 & 4$

(44)

Type: Matrix(Integer)

fricas

cost9:=trace(transpose(C9)*X9)  -- 112

$\label{eq45}112$

(45)

Type: PositiveInteger?

fricas

modiMethod(a9,b9,C9) -- check sum a = sum b

$\label{eq46}\begin{array}{@{}l} \displaystyle \left[{\hbox{\axiomType{Sol}\ } ={\left[ \begin{array}{cccc} 0 & 0 & 6 & 3 \ 0 &{12}& 3 & 0 \ {10}& 0 & 0 & 4$

(46)

Type: Record(Sol: Matrix(Integer),cost: Integer)

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