axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / MAMA.spad

)abbrev package MAMA MatrixManipulation
++ Author: Raoul Bourquin
++ Date Created: 17 November 2012
++ Date Last Updated: 1 December 2012
++ Description:
++ Some functions for manipulating (dense) matrices.
++ Supported are various kinds of slicing, splitting and stacking of 
++ matrices. The functions resemble operations often used in numerical 
++ linear algebra algorithms.

MatrixManipulation(R, Row, Col, M) : SIG == CODE where
  R : Field
  Row : FiniteLinearAggregate R
  Col : FiniteLinearAggregate R
  M : MatrixCategory(R, Row, Col)

  I ==> Integer
  PI ==> PositiveInteger
  LI ==> List I
  SI ==> Segment I
  LPI ==> List PI
  SPI ==> Segment PI

  SIG ==> with

    -- Slicing matrices

    -- How to call aRow, aColumn? Name clashed with usual row, column
    -- Package call is ugly because of many parameters of MAMA

    element : (M, PI, PI) -> M
      ++ \spad{element} returns a single element out of a matrix.
      ++ The element is put into a one by one matrix.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X element(M,2,2)

    aRow : (M, PI) -> M
      ++ \spad{aRow} returns a single row out of a matrix.
      ++ The row is put into a one by N matrix.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X aRow(M, 1)
      ++X aRow(M, 2)

    rows : (M, LPI) -> M
      ++ \spad{rows} returns several rows out of a matrix.
      ++ The rows are stacked into a matrix.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X rows(M, [1,2])
      ++X rows(M, [3,2])

    rows : (M, SPI) -> M
      ++ \spad{rows} returns several rows out of a matrix.
      ++ The rows are stacked into a matrix.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X rows(M, 2..3)

    aColumn : (M, PI) -> M
      ++ \spad{aColumn} returns a single column out of a matrix.
      ++ The column is put into a one by N matrix.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X aColumn(M, 2)

    columns : (M, LPI) -> M
      ++ \spad{columns} returns several columns out of a matrix.
      ++ The columns are stacked into a matrix.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X columns(M, [1,2])
      ++X columns(M, [3,2])

    columns : (M, SPI) -> M
      ++ \spad{columns} returns several columns out of a matrix.
      ++ The columns are stacked into a matrix.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X columns(M, 1..2)

    subMatrix : (M, LPI, LPI) -> M
      ++ \spad{subMatrix} returns several elements out of a matrix.
      ++ The elements are stacked into a submatrix.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X subMatrix(M, [1,2],[1,2])
      ++X subMatrix(M, [1,3],[1,3])

    subMatrix : (M, SPI, SPI) -> M
      ++ \spad{subMatrix} returns several elements out of a matrix.
      ++ The elements are stacked into a submatrix.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X subMatrix(M, 1..2,2..3)

    diagonalMatrix : (M, I) -> M
      ++ \spad{diagonalMatrix} returns a diagonal out of a matrix.
      ++ The diagonal is put into a matrix of same shape as the
      ++ original one. Positive integer arguments select upper
      ++ off-diagonals, negative ones lower off-diagonals.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X diagonalMatrix(M, 1)
      ++X diagonalMatrix(M, 2)
      ++X diagonalMatrix(M, -1)

    diagonalMatrix : M -> M
      ++ \spad{diagonalMatrix} returns the main diagonal out of
      ++ a matrix. The diagonal is put into a matrix of same shape
      ++ as the original one.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X diagonalMatrix(M)

    bandMatrix : (M, LI) -> M
      ++ \spad{bandMatrix} returns multiple diagonals out of a matrix.
      ++ The diagonals are put into a matrix of same shape as the
      ++ original one. Positive integer arguments select upper
      ++ off-diagonals, negative ones lower off-diagonals.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X bandMatrix(M, [-1,1])
      ++X bandMatrix(M, [-1,0,1])

    bandMatrix : (M, SI) -> M
      ++ \spad{bandMatrix} returns multiple diagonals out of a matrix.
      ++ The diagonals are put into a matrix of same shape as the
      ++ original one. Positive integer arguments select upper
      ++ off-diagonals, negative ones lower off-diagonals.
      ++
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X bandMatrix(M, -1..1)

    -- Stacking matrices

    horizConcat : (List M) -> M
      ++ \spad{horizConcat} concatenates matrices column wise.
      ++
      ++X A := matrix([[a]])
      ++X B := matrix([[b]])
      ++X C := matrix([[c]])
      ++X A12 := horizConcat([A,B,C])

    vertConcat : (List M) -> M
      ++ \spad{vertConcat} concatenates matrices row wise.
      ++
      ++X A := matrix([[a]])
      ++X B := matrix([[b]])
      ++X C := matrix([[c]])
      ++X A21 := vertConcat([A,B,C])

    blockConcat : (List List M) -> M
      ++ \spad{blockConcat} concatenates matrices row and
      ++ column wise, building a block matrix. The order
      ++ is row major as in \spad{matrix}.
      ++
      ++X A := matrix([[a]])
      ++X B := matrix([[b]])
      ++X C := matrix([[c]])
      ++X A11 := element(M, 3,3)
      ++X A12 := horizConcat([A,B,C])
      ++X A21 := vertConcat([A,B,C])
      ++X M := matrix([[a,b,c],[d,e,f],[g,h,i]])
      ++X E := blockConcat([[A11,A12],[A21,M]])
      ++X t1 := blockSplit(E, 4, [2,2])
      ++X t2 := blockConcat t1
      ++X zero?(E-t2)
      ++X t3 := blockSplit(E, [1,2,1], [2,2])
      ++X t4 := blockConcat t3
      ++X zero?(E-t4)

    -- Splitting matrices

    vertSplit : (M, PI) -> List M
      ++ \spad{vertSplit} splits a matrix into multiple
      ++ submatrices row wise.
      ++
      ++X E := matrix([[i,a,b,c],[a,a,b,c],[b,d,e,f],[c,g,h,i]])
      ++X t1:= vertSplit(E, 2)

    vertSplit : (M, LPI) -> List M
      ++ \spad{vertSplit} splits a matrix into multiple
      ++ submatrices row wise.
      ++
      ++X E := matrix([[i,a,b,c],[a,a,b,c],[b,d,e,f],[c,g,h,i]])
      ++X t1:= vertSplit(E, [1,2,1])

    horizSplit : (M, PI) -> List M
      ++ \spad{horizSplit} splits a matrix into multiple
      ++ submatrices column wise.
      ++
      ++X E := matrix([[i,a,b,c],[a,a,b,c],[b,d,e,f],[c,g,h,i]])
      ++X t1:= horizSplit(E, 2)

    horizSplit : (M, LPI) -> List M
      ++ \spad{horizSplit} splits a matrix into multiple
      ++ submatrices column wise.
      ++
      ++X E := matrix([[i,a,b,c],[a,a,b,c],[b,d,e,f],[c,g,h,i]])
      ++X t1:= horizSplit(E, [2,2])
      ++X t2:= horizSplit(E, [1,2,1])

    blockSplit : (M, PI, PI) -> List List M
      ++ \spad{blockSplit} splits a matrix into multiple
      ++ submatrices row and column wise, dividing
      ++ a matrix into blocks.
      ++
      ++X E := matrix([[i,a,b,c],[a,a,b,c],[b,d,e,f],[c,g,h,i]])
      ++X t1:= blockSplit(E,2,2)

    blockSplit : (M, LPI, PI) -> List List M
      ++ \spad{blockSplit} splits a matrix into multiple
      ++ submatrices row and column wise, dividing
      ++ a matrix into blocks.
      ++
      ++X E := matrix([[i,a,b,c],[a,a,b,c],[b,d,e,f],[c,g,h,i]])
      ++X t1:= blockSplit(E, [2,1,1], 2)

    blockSplit : (M, PI, LPI) -> List List M
      ++ \spad{blockSplit} splits a matrix into multiple
      ++ submatrices row and column wise, dividing
      ++ a matrix into blocks.
      ++
      ++X E := matrix([[i,a,b,c],[a,a,b,c],[b,d,e,f],[c,g,h,i]])
      ++X t1:= blockSplit(E, 4, [2,2])

    blockSplit : (M, LPI, LPI) -> List List M
      ++ \spad{blockSplit} splits a matrix into multiple
      ++ submatrices row and column wise, dividing
      ++ a matrix into blocks.
      ++
      ++X E := matrix([[i,a,b,c],[a,a,b,c],[b,d,e,f],[c,g,h,i]])
      ++X t1:= blockSplit(E, [1,2,1], [2,2])

  CODE ==> add

    minr ==> minRowIndex
    maxr ==> maxRowIndex
    minc ==> minColIndex
    maxc ==> maxColIndex

    -- Custom function to expand Segment(PositiveInteger) into 
    -- List(PositiveInteger). This operation is not supported by the 
    -- overly restrictive library implementation.
    expand(spi : SPI) : LPI ==
        lr := empty()$LPI
        l : PI := lo spi
        h : PI := hi spi
        inc : I := incr spi
        zero? inc => error "Cannot expand a segment with an increment of zero"
        if inc > 0 then
          while l <= h repeat
            lr := concat(l, lr)
            l := (l + inc) pretend PI
        else
          while l >= h repeat
            lr := concat(l, lr)
            l := (l + inc) pretend PI
        reverse! lr

    element(A, r, c) ==
      matrix([[A(r,c)]])

    aRow(A:M, r:PI) : M ==
      subMatrix(A, r, r, minc A, maxc A)

    rows(A:M, lst:LPI) : M ==
      ls := [aRow(A, r) for r in lst]
      reduce(vertConcat, ls)

    rows(A:M, si:SPI) : M ==
      rows(A, expand(si))

    aColumn(A:M, c:PI) : M ==
      subMatrix(A, minr A, maxr A, c, c)

    columns(A:M, lst:LPI) : M ==
      ls := [aColumn(A,c) for c in lst]
      reduce(horizConcat, ls)

    columns(A:M, si:SPI) : M ==
      columns(A, expand(si))

    diagonalMatrix(A, n) ==
      nr := nrows(A)
      nc := ncols(A)
      n > (nc-1) => error "requested diagonal out of range"
      n < 0 and abs(n) > (nr-1) => error "requested diagonal out of range"
      B := zero(nr,nc)
      if n >= 0 then
        dl := min(nc-n, nr)
        sr := minr(A)
        sc := minc(A) + n
      else
        dl := min(nc, nr-abs(n))
        sr := minr(A) + abs(n)
        sc := minc(A)
      for i in 0..(dl-1) repeat
        qsetelt!(B, sr+i, sc+i, A(sr+i, sc+i))
      B

    diagonalMatrix(A) ==
      diagonalMatrix(A, 0)

    bandMatrix(A:M, ln:LI) : M ==
      -- Really inefficient
      reduce("+", [diagonalMatrix(A,d) for d in ln])

    bandMatrix(A:M, si:SI) : M ==
      bandMatrix(A, expand(si))

    subMatrix(A:M, lr:LPI, lc:LPI) : M ==
      -- Really inefficient
      lle := [[ element(A,r,c) for c in lc] for r in lr]
      blockConcat(lle)

    subMatrix(A:M, sr:SPI, sc:SPI) : M ==
      subMatrix(A, low sr, high sr, low sc, high sc)

    -- Stack matrices

    horizConcat(LA) ==
      reduce(horizConcat, LA)

    vertConcat(LA) ==
      reduce(vertConcat, LA)

    blockConcat(LLA: List List M) : M ==
      reduce(vertConcat, [reduce(horizConcat, LA) for LA in LLA])

    -- Split matrices

    vertSplit(A:M, r:PI) : List M ==
      dr := nrows(A) exquo r
      dr case "failed" => error "split does not result in an equal division"
      mir := minr A
      mic := minc A
      mac := maxc A
      [ subMatrix(A, mir+i*dr, mir+(i+1)*dr-1, mic, mac) for i in 0..(r-1) ]

    vertSplit(A:M, lr:LPI) : List M ==
      reduce("+", lr) ~= nrows(A) => _
          error "split does not result in proper partition"
      l : List PI := cons(1, scan(_+, lr, 1$PI)$ListFunctions2(PI,PI))
      mir := minr(A) -1   -- additional shift because l starts at 1
      mic := minc A
      mac := maxc A
      result := _
        [ subMatrix(A, mir+l(i-1), mir+l(i)-1, mic, mac) for i in 2..#l ]

    horizSplit(A:M, c:PI) : List M ==
      dc := ncols(A) exquo c
      dc case "failed" => error "split does not result in an equal division"
      mir := minr A
      mar := maxr A
      mic := minc A
      [ subMatrix(A, mir, mar, mic+i*dc, mic+(i+1)*dc-1) for i in 0..(c-1) ]

    horizSplit(A:M, lc:LPI) : List M ==
      reduce("+", lc) ~= ncols(A) => _
          error "split does not result in proper partition"
      l : List PI := cons(1, scan(_+, lc, 1$PI)$ListFunctions2(PI,PI))
      mir := minr A
      mar := maxr A
      mic := minc(A) -1   -- additional shift because l starts at 1
      result := _
         [ subMatrix(A, mir, mar, mic+l(i-1), mic+l(i)-1) for i in 2..#l ]

    blockSplit(A:M, nr:PI, nc:PI) : List List M ==
      [ horizSplit(X, nc) for X in vertSplit(A, nr) ]

    blockSplit(A:M, lr:LPI, nc:PI) : List List M ==
      [ horizSplit(X, nc) for X in vertSplit(A, lr) ]

    blockSplit(A:M, nr:PI, lc:LPI) : List List M ==
      [ horizSplit(X, lc) for X in vertSplit(A, nr) ]

    blockSplit(A:M, lr:LPI, lc:LPI) : List List M ==
      [ horizSplit(X, lc) for X in vertSplit(A, lr) ]