|
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]
spad 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.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
"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 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/DOPTTest flavours
fricas --)co dopt
- supply, demand
a:=vector [4,7,4]$Vector INT
Type: Vector(Integer)
fricas b:=vector [2,4,6,3]$Vector INT
Type: Vector(Integer)
fricas -- cost matrix
C:= matrix [[2,3,4,1],[5,4,2,3],[4,2,8,6]]
Type: Matrix(Integer)
fricas X:=computeSolutionMatrix(3,4,a,b,C)
Type: Matrix(Integer)
fricas cost:=trace(transpose(C)*X) -- 29
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]]
Type: Matrix(Integer)
fricas b2:= vector [ 30, 50, 20, 40, 30, 11]
Type: Vector(PositiveInteger ?)
fricas a2:= vector [ 50, 40, 60, 31]
Type: Vector(PositiveInteger ?)
fricas X2:=computeSolutionMatrix(4,6,a2,b2,C2)
Type: Matrix(Integer)
fricas cost2:=trace(transpose(C2)*X2) -- 330
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]]
Type: Matrix(Integer)
fricas b3:= vector [ 20, 40, 30, 10, 50, 25]
Type: Vector(PositiveInteger ?)
fricas a3:= vector [ 30, 50, 75, 20]
Type: Vector(PositiveInteger ?)
fricas X3:=computeSolutionMatrix(4,6,a3,b3,C3)
Type: Matrix(Integer)
fricas cost3:=trace(transpose(C3)*X3) -- 430
fricas -- 4 ??? sum(a)~sum(b)
C4:= matrix [[ 5, 3, 6, 2],[ 4, 7, 9, 1],[ 3, 4, 7, 5]]
Type: Matrix(Integer)
fricas b4 := vector [ 16, 18, 30, 25]
Type: Vector(PositiveInteger ?)
fricas a4 := vector [ 19, 37, 34]
Type: Vector(PositiveInteger ?)
fricas X4:=computeSolutionMatrix(3,4,a4,b4,C4)
Type: Matrix(Integer)
fricas cost4:=trace(transpose(C4)*X4) -- 348
fricas -- 5 WP TT
C5:= matrix [[6,4,3],[3,2,2]]
Type: Matrix(Integer)
fricas a5:= vector [16,14]
Type: Vector(PositiveInteger ?)
fricas b5:= vector [15,5,10]
Type: Vector(PositiveInteger ?)
fricas X5:=computeSolutionMatrix(2,3,a5,b5,C5)
Type: Matrix(Integer)
fricas cost5:=trace(transpose(C5)*X5) -- 98
fricas -- 6 WP ZykM
C6:= matrix [[7,4,3],[5,5,6]]
Type: Matrix(Integer)
fricas a6:= vector [12,8]
Type: Vector(PositiveInteger ?)
fricas b6:= vector [4,10,6]
Type: Vector(PositiveInteger ?)
fricas X6:=computeSolutionMatrix(2,3,a6,b6,C6)
Type: Matrix(Integer)
fricas cost6:=trace(transpose(C6)*X6) -- 82
fricas --- 7
C7:= matrix [[140,130,170,150,140],[150,160,140,180,130],[160,120,180,170,190]]
Type: Matrix(Integer)
fricas a7:= vector [900,800,400]
Type: Vector(PositiveInteger ?)
fricas b7:= vector [300,500,400,600,300]
Type: Vector(PositiveInteger ?)
fricas X7:=computeSolutionMatrix(3,5,a7,b7,C7)
Type: Matrix(Integer)
fricas cost7:=trace(transpose(C7)*X7) -- 289000
fricas --- 8
C8:= matrix [[3,2],[1,5],[5,4]]
Type: Matrix(Integer)
fricas a8:= vector [45,60,5]
Type: Vector(PositiveInteger ?)
fricas b8:= vector [50,60]
Type: Vector(PositiveInteger ?)
fricas X8:=computeSolutionMatrix(3,2,a8,b8,C8)
Type: Matrix(Integer)
fricas cost8:=trace(transpose(C8)*X8) -- 210
fricas --- 9
C9:= matrix [[5,5,2,7],[8,3,1,8],[2,4,4,5]]
Type: Matrix(Integer)
fricas a9:= vector [9,15,14]
Type: Vector(PositiveInteger ?)
fricas b9:= vector [10,12,9,7]
Type: Vector(PositiveInteger ?)
fricas X9:=computeSolutionMatrix(3,4,a9,b9,C9)
Type: Matrix(Integer)
fricas cost9:=trace(transpose(C9)*X9) -- 112
fricas modiMethod(a9,b9,C9) -- check sum a = sum b
Type: Record(Sol: Matrix(Integer),cost: Integer)
Note: "null" must have been removed recently ??
WA: null ==> empty?
|
![\label{eq1}\left[ 4, \: 7, \: 4 \right]](images/2293666038040784793-16.0px.png "\label{eq1}\left[ 4, \: 7, \: 4 \right]")
![\label{eq2}\left[ 2, \: 4, \: 6, \: 3 \right]](images/1288975305893543674-16.0px.png "\label{eq2}\left[ 2, \: 4, \: 6, \: 3 \right]")
![\label{eq3}\left[
\begin{array}{cccc}
2 & 3 & 4 & 1
\
5 & 4 & 2 & 3
\
4 & 2 & 8 & 6](images/307782075778182554-16.0px.png "\label{eq3}\left[
\begin{array}{cccc}
2 & 3 & 4 & 1
\
5 & 4 & 2 & 3
\
4 & 2 & 8 & 6")
![\label{eq4}\left[
\begin{array}{cccc}
2 & 0 & 0 & 2
\
0 & 0 & 6 & 1
\
0 & 4 & 0 & 0](images/922101987126892618-16.0px.png "\label{eq4}\left[
\begin{array}{cccc}
2 & 0 & 0 & 2
\
0 & 0 & 6 & 1
\
0 & 4 & 0 & 0")
![\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](images/8486860699344650811-16.0px.png "\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")
![\label{eq7}\left[{30}, \:{50}, \:{20}, \:{40}, \:{30}, \:{11}\right]](images/496489830757336108-16.0px.png "\label{eq7}\left[{30}, \:{50}, \:{20}, \:{40}, \:{30}, \:{11}\right]")
![\label{eq8}\left[{50}, \:{40}, \:{60}, \:{31}\right]](images/8824075216539292821-16.0px.png "\label{eq8}\left[{50}, \:{40}, \:{60}, \:{31}\right]")
![\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](images/4265692357471623946-16.0px.png "\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")
![\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](images/7818496027194715769-16.0px.png "\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")
![\label{eq12}\left[{20}, \:{40}, \:{30}, \:{10}, \:{50}, \:{25}\right]](images/2632458141999258084-16.0px.png "\label{eq12}\left[{20}, \:{40}, \:{30}, \:{10}, \:{50}, \:{25}\right]")
![\label{eq13}\left[{30}, \:{50}, \:{75}, \:{20}\right]](images/4330086038618137733-16.0px.png "\label{eq13}\left[{30}, \:{50}, \:{75}, \:{20}\right]")
![\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](images/2661397530347406787-16.0px.png "\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")
![\label{eq16}\left[
\begin{array}{cccc}
5 & 3 & 6 & 2
\
4 & 7 & 9 & 1
\
3 & 4 & 7 & 5](images/8068846246271712093-16.0px.png "\label{eq16}\left[
\begin{array}{cccc}
5 & 3 & 6 & 2
\
4 & 7 & 9 & 1
\
3 & 4 & 7 & 5")
![\label{eq17}\left[{16}, \:{18}, \:{30}, \:{25}\right]](images/1875562298789017105-16.0px.png "\label{eq17}\left[{16}, \:{18}, \:{30}, \:{25}\right]")
![\label{eq18}\left[{19}, \:{37}, \:{34}\right]](images/3217710673625623977-16.0px.png "\label{eq18}\left[{19}, \:{37}, \:{34}\right]")
![\label{eq19}\left[
\begin{array}{cccc}
0 &{18}& 1 & 0
\
{12}& 0 & 0 &{25}
\
4 & 0 &{29}& 0](images/5036562955384049334-16.0px.png "\label{eq19}\left[
\begin{array}{cccc}
0 &{18}& 1 & 0
\
{12}& 0 & 0 &{25}
\
4 & 0 &{29}& 0")
![\label{eq21}\left[
\begin{array}{ccc}
6 & 4 & 3
\
3 & 2 & 2](images/8684554170182937113-16.0px.png "\label{eq21}\left[
\begin{array}{ccc}
6 & 4 & 3
\
3 & 2 & 2")
![\label{eq22}\left[{16}, \:{14}\right]](images/8751500832607147413-16.0px.png "\label{eq22}\left[{16}, \:{14}\right]")
![\label{eq23}\left[{15}, \: 5, \:{10}\right]](images/2097084593947777502-16.0px.png "\label{eq23}\left[{15}, \: 5, \:{10}\right]")
![\label{eq24}\left[
\begin{array}{ccc}
1 & 5 &{10}
\
{14}& 0 & 0](images/5270394019479351289-16.0px.png "\label{eq24}\left[
\begin{array}{ccc}
1 & 5 &{10}
\
{14}& 0 & 0")
![\label{eq26}\left[
\begin{array}{ccc}
7 & 4 & 3
\
5 & 5 & 6](images/2774735784986788120-16.0px.png "\label{eq26}\left[
\begin{array}{ccc}
7 & 4 & 3
\
5 & 5 & 6")
![\label{eq27}\left[{12}, \: 8 \right]](images/3905627965177140290-16.0px.png "\label{eq27}\left[{12}, \: 8 \right]")
![\label{eq28}\left[ 4, \:{10}, \: 6 \right]](images/5578750935806267611-16.0px.png "\label{eq28}\left[ 4, \:{10}, \: 6 \right]")
![\label{eq29}\left[
\begin{array}{ccc}
0 & 6 & 6
\
4 & 4 & 0](images/1184337983058265501-16.0px.png "\label{eq29}\left[
\begin{array}{ccc}
0 & 6 & 6
\
4 & 4 & 0")
![\label{eq31}\left[
\begin{array}{ccccc}
{140}&{130}&{170}&{150}&{140}
\
{150}&{160}&{140}&{180}&{130}
\
{160}&{120}&{180}&{170}&{190}](images/3575101501529275667-16.0px.png "\label{eq31}\left[
\begin{array}{ccccc}
{140}&{130}&{170}&{150}&{140}
\
{150}&{160}&{140}&{180}&{130}
\
{160}&{120}&{180}&{170}&{190}")
![\label{eq32}\left[{900}, \:{800}, \:{400}\right]](images/2762589329817187818-16.0px.png "\label{eq32}\left[{900}, \:{800}, \:{400}\right]")
![\label{eq33}\left[{300}, \:{500}, \:{400}, \:{600}, \:{300}\right]](images/4085208705409753303-16.0px.png "\label{eq33}\left[{300}, \:{500}, \:{400}, \:{600}, \:{300}\right]")
![\label{eq34}\left[
\begin{array}{ccccc}
{200}&{100}& 0 &{600}& 0
\
{100}& 0 &{400}& 0 &{300}
\
0 &{400}& 0 & 0 & 0](images/6740782615661095736-16.0px.png "\label{eq34}\left[
\begin{array}{ccccc}
{200}&{100}& 0 &{600}& 0
\
{100}& 0 &{400}& 0 &{300}
\
0 &{400}& 0 & 0 & 0")
![\label{eq36}\left[
\begin{array}{cc}
3 & 2
\
1 & 5
\
5 & 4](images/1321132779521996323-16.0px.png "\label{eq36}\left[
\begin{array}{cc}
3 & 2
\
1 & 5
\
5 & 4")
![\label{eq37}\left[{45}, \:{60}, \: 5 \right]](images/7450196921999374166-16.0px.png "\label{eq37}\left[{45}, \:{60}, \: 5 \right]")
![\label{eq38}\left[{50}, \:{60}\right]](images/4262996343123119087-16.0px.png "\label{eq38}\left[{50}, \:{60}\right]")
![\label{eq39}\left[
\begin{array}{cc}
0 &{45}
\
{50}&{10}
\
0 & 5](images/1572605817629471512-16.0px.png "\label{eq39}\left[
\begin{array}{cc}
0 &{45}
\
{50}&{10}
\
0 & 5")
![\label{eq41}\left[
\begin{array}{cccc}
5 & 5 & 2 & 7
\
8 & 3 & 1 & 8
\
2 & 4 & 4 & 5](images/6976500074622694705-16.0px.png "\label{eq41}\left[
\begin{array}{cccc}
5 & 5 & 2 & 7
\
8 & 3 & 1 & 8
\
2 & 4 & 4 & 5")
![\label{eq42}\left[ 9, \:{15}, \:{14}\right]](images/4724158966728207473-16.0px.png "\label{eq42}\left[ 9, \:{15}, \:{14}\right]")
![\label{eq43}\left[{10}, \:{12}, \: 9, \: 7 \right]](images/6821201727280907159-16.0px.png "\label{eq43}\left[{10}, \:{12}, \: 9, \: 7 \right]")
![\label{eq44}\left[
\begin{array}{cccc}
0 & 0 & 6 & 3
\
0 &{12}& 3 & 0
\
{10}& 0 & 0 & 4](images/320298404747872023-16.0px.png "\label{eq44}\left[
\begin{array}{cccc}
0 & 0 & 6 & 3
\
0 &{12}& 3 & 0
\
{10}& 0 & 0 & 4")
![\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](images/5355187779144050539-16.0px.png "\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")