http://mathling.com/core/array  library module

http://mathling.com/core/array


2D array of numbers: similar functionality to matrix.xqy, but with
sequences of doubles. Similar to the relationship between point and vector.
This supports more matrix operations, such as transposition and inverse.
This is a row major ordering: note that the get() function, like those
for point matrices of various kinds puts the row first e.g. get(y,x)
which is, I admit a bit confusing.

Copyright© Mary Holstege 2021-2023
CC-BY (https://creativecommons.org/licenses/by/4.0/)

December 2021
Status: Stable

Imports

http://mathling.com/core/utilities
import module namespace util="http://mathling.com/core/utilities"
       at "../core/utilities.xqy"
http://mathling.com/core/config
import module namespace config="http://mathling.com/core/config"
       at "../core/config.xqy"
http://mathling.com/core/errors
import module namespace errors="http://mathling.com/core/errors"
       at "../core/errors.xqy"

Variables

Variable: $RESERVED as 

Functions

Function: array
declare function array($rows as xs:integer, $columns as xs:integer, $data as xs:double*) as map(*)


array()
Construct a new matrix

Params
  • rows as xs:integer number of rows in matrix
  • columns as xs:integer number of columns in matrix
  • data as xs:double* data values (must be $rows X $columns)
Returns
  • map(*): matrix
declare function this:array(
  $rows as xs:integer,
  $columns as xs:integer,
  $data as xs:double*
) as map(*)
{
  if ($rows * $columns ne count($data))
  then errors:error("UTIL-BADARRAY", ($rows, $columns, count($data)))
  else (
    map {
      "kind": "array",
      "rows": $rows,
      "columns": $columns,
      "data": $data
    }
  )
}
Errors

UTIL-BADARRAY if number of data values doesn't match size

Function: rows
declare function rows($array as map(*)) as xs:integer


rows()
Accessor for the number of rows in the array

Params
  • array as map(*) the array
Returns
  • xs:integer: number of rows in array
declare function this:rows($array as map(*)) as xs:integer
{
  $array("rows")
}

Function: columns
declare function columns($array as map(*)) as xs:integer


columns()
Accessor for the number of rows in the array

Params
  • array as map(*) the array
Returns
  • xs:integer: number of rows in array
declare function this:columns($array as map(*)) as xs:integer
{
  $array("columns")
}

Function: data
declare function data($array as map(*)) as xs:double*


data()
Accessor for the data in the array

Params
  • array as map(*) the array
Returns
  • xs:double*: sequence of values (row major order)
declare function this:data($array as map(*)) as xs:double*
{
  $array("data")
}

Function: data
declare function data($array as map(*), $data as xs:double*) as map(*)


data()
Settor for the data in the array

Params
  • array as map(*) the array
  • data as xs:double* the new data
Returns
  • map(*): new matrix
declare function this:data($array as map(*), $data as xs:double*) as map(*)
{
  if (this:rows($array) * this:columns($array) ne count($data))
  then errors:error("UTIL-BADARRAY", (this:rows($array), this:columns($array), count($data)))
  else (
    $array=>map:put("data", $data)
  )
}
Errors

UTIL-BADARRAY if number of data values doesn't match size

Function: get
declare function get($array as map(*), $row as xs:integer, $column as xs:integer) as xs:double


get()
Get a specific data value

Params
  • array as map(*): the array
  • row as xs:integer: the row number
  • column as xs:integer: the column number
Returns
  • xs:double: data value array[row, column]
declare function this:get(
  $array as map(*),
  $row as xs:integer,
  $column as xs:integer
) as xs:double
{
  this:data($array)[($row - 1)*$array("columns") + $column]
}

Function: row
declare function row($array as map(*), $row as xs:integer) as xs:double*


row()
Return all the values from the given row

Params
  • array as map(*): the array
  • row as xs:integer: the row number
Returns
  • xs:double*: values in row in column order
declare function this:row(
  $array as map(*),
  $row as xs:integer
) as xs:double*
{
  if ($row lt 1 or $row gt $array("rows")) then () else
  let $data := $array("data")
  let $num-columns := $array("columns")
  let $ixs :=
    for $col in 1 to $num-columns
    return ($row - 1)*$num-columns + $col
  for $ix in $ixs
  return $data[$ix]
}

Function: column
declare function column($array as map(*), $column as xs:integer) as xs:double*


column()
Return all the values from the given column

Params
  • array as map(*): the array
  • column as xs:integer: the column number
Returns
  • xs:double*: values in column in row order
declare function this:column(
  $array as map(*),
  $column as xs:integer
) as xs:double*
{
  if ($column lt 1 or $column gt $array("columns")) then () else
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  let $ixs :=
    for $row in 1 to $num-rows
    return ($row - 1)*$num-columns + $column
  for $ix in $ixs
  return $data[$ix]
}

Function: diagonal
declare function diagonal($array as map(*)) as xs:double*


diagonal()
Return all the values from the diagonal
Raises an error if the array is not square

Params
  • array as map(*): the array
Returns
  • xs:double*: values along main diagonal
declare function this:diagonal($array as map(*)) as xs:double*
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    if ($num-rows != $num-columns) then (
      errors:error("UTIL-MATRIXNOTSQUARE", ($num-rows,$num-columns))
    ) else (
      let $ixs := for $i in 1 to $num-rows return ($i - 1)*$num-columns + $i
      for $ix in $ixs
      return $data[$ix]
    )
  )
}
Errors

UTIL-MATRIXNOTSQUARE if matrix is not square

Function: trace
declare function trace($array as map(*)) as xs:double


trace()
Compute the trace of the array.
Raises an error if the array is not square

Params
  • array as map(*): the array
Returns
  • xs:double: trace of the array
declare function this:trace($array as map(*)) as xs:double
{
  sum(this:diagonal($array))
}
Errors

UTIL-MATRIXNOTSQUARE if matrix is not square

Function: determinant
declare function determinant($array as map(*)) as xs:double


determinant()
Compute the determinant of the array.
Raises an error if the array is not square

Params
  • array as map(*): the array
Returns
  • xs:double: determinant
declare function this:determinant($array as map(*)) as xs:double
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    if ($num-rows != $num-columns) then (
      errors:error("UTIL-MATRIXNOTSQUARE", ($num-rows,$num-columns))
    ) else (
      switch($num-rows)
      case 1 return $data
      case 2 return $data[1]*$data[4] - $data[2]*$data[3]
      case 3 return (
         $data[1]*$data[5]*$data[9] + $data[2]*$data[6]*$data[7] + $data[3]*$data[4]*$data[8]
        -$data[1]*$data[6]*$data[8] - $data[2]*$data[4]*$data[9] - $data[3]*$data[5]*$data[7]
      )
      default return (
        sum(
          for $i in 1 to $num-columns return (
            math:pow(-1, 1 + $i) *
            this:determinant(
              this:array(
                $num-rows - 1,
                $num-columns - 1,
                for $r in 1 to $num-rows
                for $c in 1 to $num-columns
                where $r != 1 and $c != $i
                return $data[($r - 1)*$num-columns + $c]
              )
            )
          )
        )
      )
    )
  )
}
Errors

UTIL-MATRIXNOTSQUARE if matrix is not square

Function: transpose
declare function transpose($array as map(*)) as map(*)


transpose()
Transpose the array.

Params
  • array as map(*): the array
Returns
  • map(*): the array, transposed
declare function this:transpose($array as map(*)) as map(*)
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    this:array(
      $num-columns,
      $num-rows,
      for $c in 1 to $num-columns
      for $r in 1 to $num-rows
      return $data[($r - 1)*$num-columns + $c]
    )
  )
}

Function: cofactor-matrix
declare function cofactor-matrix($array as map(*)) as map(*)


cofactor-matrix()
Compute the cofactor matrix for the array. Raises an error if it is not square.

Params
  • array as map(*): the array
Returns
  • map(*): matrix of cofactors
declare function this:cofactor-matrix($array as map(*)) as map(*)
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    if ($num-rows != $num-columns) then (
      errors:error("UTIL-MATRIXNOTSQUARE", ($num-rows,$num-columns))
    ) else if ($num-rows = 1) then (
      errors:error("UTIL-CANNOT", ("cofactor-matrix", $data))
    ) else (
      this:array(
        $num-rows,
        $num-columns,
        for $j in 1 to $num-rows
        for $i in 1 to $num-columns
        return (
          math:pow(-1, $i + $j) *
          this:determinant(
            this:array(
              $num-rows - 1,
              $num-columns - 1,
              for $r in 1 to $num-rows
              for $c in 1 to $num-columns
              where $r != $j and $c != $i
              return $data[($r - 1)*$num-columns + $c]
            )
          )
        )
      )
    )
  )
}
Errors

UTIL-MATRIXNOTSQUARE if matrix is not square

Errors

UTIL-CANNOT if matrix is only 1x1

Function: adjoint
declare function adjoint($array as map(*)) as map(*)


adjoint()
Compute the adjoint of the array. Raises an error if it is not square.

Params
  • array as map(*): the array
Returns
  • map(*): adjoint matrix
declare function this:adjoint($array as map(*)) as map(*)
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    if ($num-rows != $num-columns) then (
      errors:error("UTIL-MATRIXNOTSQUARE", ($num-rows,$num-columns))
    ) else (
      this:transpose(this:cofactor-matrix($array))
    )
  )
}
Errors

UTIL-MATRIXNOTSQUARE if matrix is not square

Function: inverse
declare function inverse($array as map(*)) as map(*)


inverse()
Compute the inverse of the array. Raises an error if it is not square.

Params
  • array as map(*): the array
Returns
  • map(*): array inverse
declare function this:inverse($array as map(*)) as map(*)
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    if ($num-rows != $num-columns) then (
      errors:error("UTIL-MATRIXNOTSQUARE", ($num-rows,$num-columns))
    ) else (
      let $det := this:determinant($array)
      let $adj := this:adjoint($array)
      return this:map($adj, function ($c as xs:double) as xs:double {$c div $det})
    )
  )
}
Errors

UTIL-MATRIXNOTSQUARE if matrix is not square

Function: times
declare function times($array as map(*), $k as xs:double) as map(*)


times()
Scalar multiplication of an array by a constant.

Params
  • array as map(*): the array
  • k as xs:double: multiplier
Returns
  • map(*): k*array
declare function this:times($array as map(*), $k as xs:double) as map(*)
{
  this:map($array, function ($c as xs:double) { $c*$k })
}

Function: multiply
declare function multiply($a1 as map(*), $a2 as map(*)) as map(*)


multiply()
Matrix multiplication.

Params
  • a1 as map(*): one array
  • a2 as map(*): another array
Returns
  • map(*): a1*a2
declare function this:multiply($a1 as map(*), $a2 as map(*)) as map(*)
{
  if (this:columns($a1) != this:rows($a2)) then (
    errors:error("UTIL-MATRIXNOTCOMPAT", ("multiply", this:rows($a1), this:columns($a1), this:rows($a2), this:columns($a2)))
  ) else (
    let $rows1 := this:rows($a1)
    let $columns1 := this:columns($a1)
    let $data1 := this:data($a1)
    let $rows2 := this:rows($a2)
    let $columns2 := this:columns($a2)
    let $data2 := this:data($a2)
    return (
      this:array($rows1, $columns2, 
        for $r in 1 to $rows1
        for $c in 1 to $columns2
        let $row := this:row($a1, $r)
        let $col := this:column($a2, $c)
        return sum(for $d at $i in $row return $d * $col[$i])
      )
    )
  )
}
Errors

UTIL-MATRIXNOTCOMPAT if the array sizes are incompatible

Function: add
declare function add($a1 as map(*), $a2 as map(*)) as map(*)


add()
Matrix addition

Params
  • a1 as map(*): one array
  • a2 as map(*): another array
Returns
  • map(*): a1+a2
declare function this:add($a1 as map(*), $a2 as map(*)) as map(*)
{
  if (this:columns($a1) != this:columns($a2) or this:rows($a1) != this:rows($a2)) then (
    errors:error("UTIL-MATRIXNOTCOMPAT", ("add", this:rows($a1), this:columns($a1), this:rows($a2), this:columns($a2)))
  ) else (
    let $rows1 := this:rows($a1)
    let $columns1 := this:columns($a1)
    let $data1 := this:data($a1)
    let $rows2 := this:rows($a2)
    let $columns2 := this:columns($a2)
    let $data2 := this:data($a2)
    return (
      this:array($rows1, $columns1, 
        for $d at $i in $data1 return $d + $data2[$i]
      )
    )
  )
}
Errors

UTIL-MATRIXNOTCOMPAT if the array sizes are incompatible

Function: subtract
declare function subtract($a1 as map(*), $a2 as map(*)) as map(*)


subtract()
Matrix subtraction

Params
  • a1 as map(*): one array
  • a2 as map(*): another array
Returns
  • map(*): a1-a2
declare function this:subtract($a1 as map(*), $a2 as map(*)) as map(*)
{
  if (this:columns($a1) != this:columns($a2) or this:rows($a1) != this:rows($a2)) then (
    errors:error("UTIL-MATRIXNOTCOMPAT", ("subtract", this:rows($a1), this:columns($a1), this:rows($a2), this:columns($a2)))
  ) else (
    let $rows1 := this:rows($a1)
    let $columns1 := this:columns($a1)
    let $data1 := this:data($a1)
    let $rows2 := this:rows($a2)
    let $columns2 := this:columns($a2)
    let $data2 := this:data($a2)
    return (
      this:array($rows1, $columns1, 
        for $d at $i in $data1 return $d - $data2[$i]
      )
    )
  )
}
Errors

UTIL-MATRIXNOTCOMPAT if the array sizes are incompatible

Function: map
declare function map($array as map(*), $f as function(xs:double) as xs:double) as map(*)


map()
Run a function over every item in the array, return the new array, preserving properties.

Params
  • array as map(*): the array
  • f as function(xs:double)asxs:double: function to apply to all the data in the array
Returns
  • map(*): new array
declare function this:map($array as map(*), $f as function(xs:double) as xs:double) as map(*)
{
  util:merge-into(
    $array,
    map {"data": for-each(this:data($array), $f)}
  )
}

Function: snap
declare function snap($array as map(*)) as map(*)


snap()
Snap every value in the array to an integer. Returns an array with integer data values.
Preserves properties.

Params
  • array as map(*): the array
Returns
  • map(*): new array
declare function this:snap($array as map(*)) as map(*)
{
  util:merge-into(
    $array,
    map {
      "data":
        for-each(this:data($array),
          function ($c as xs:double) as xs:integer {
            fn:round($c) cast as xs:integer
          }
        )
    }
  )
}

Function: floor
declare function floor($array as map(*)) as map(*)


floor()
Snap every value in the array to the floor of that value, as an integer.
Returns an array with integer data values. Preserves properties.

Params
  • array as map(*): the array
Returns
  • map(*): new array
declare function this:floor($array as map(*)) as map(*)
{
  util:merge-into(
    $array,
    map {
      "data": 
        for-each(this:data($array), 
          function($c as xs:double) as xs:integer {
            fn:floor($c) cast as xs:integer
          }
        )
    }
  )
}

Function: decimal
declare function decimal($array as map(*), $digits as xs:integer) as map(*)


decimal()
Truncate the number of digits in every value in the array to the given number of digits,
returning the new array. Preserves properties.

Params
  • array as map(*): the array
  • digits as xs:integer: number of digits to preserve
Returns
  • map(*): new array
declare function this:decimal($array as map(*), $digits as xs:integer) as map(*)
{
  util:merge-into(
    $array,
    map {
      "data":
        for-each(this:data($array),
          function($c as xs:double) as xs:double {
            util:decimal($c, $digits)
          }
        )
    }
  )
}

Function: quote
declare function quote($array as map(*)) as xs:string


quote()
Return a string representation of the array, suitable for debugging.

Params
  • array as map(*): the array
Returns
  • xs:string: string
declare function this:quote($array as map(*)) as xs:string
{
  string-join((
    "[",
    for $row in 1 to this:rows($array)
    return string-join(this:row($array,$row),","),
    "]"
  ),"
")
}

Function: same
declare function same($this as map(*), $other as map(*)) as xs:boolean


same()
Equality test for two arrays. Arrays are the same if they have the same dimensions and data.
Auxiliary properties do not factor into this test.

Params
  • this as map(*): one array
  • other as map(*)
Returns
  • xs:boolean: true if the arrays have the same data and sizes
declare function this:same($this as map(*), $other as map(*)) as xs:boolean
{
  this:rows($this)=this:rows($other) and
  this:columns($this)=this:columns($other) and
  deep-equal(this:data($this), this:data($other))
}

Function: horizontal-blur
declare function horizontal-blur($source as map(*), $radius as xs:integer) as map(*)


horizontal-blur()
Quick horizontal blurring.

Params
  • source as map(*) source array
  • radius as xs:integer radius of blurring (i.e number of pixels distant)
Returns
  • map(*): array with values subject to horizontal blurring
declare function this:horizontal-blur (
  $source as map(*),
  $radius as xs:integer
) as map(*)
{
  let $rows := this:rows($source)
  let $columns := this:columns($source)
  let $data := this:data($source)

  let $dest-data := (
    for $y in 1 to this:rows($source)
    let $total1 :=
      sum(for $kx in -$radius to $radius return $data[($y - 1)*$columns + $kx])

    let $totals :=
      fold-left(2 to $columns, $total1,
        function($totals as xs:double*, $x as xs:integer) {
          $totals,
          $totals[last()]
            - ($data[($y - 1)*$columns + $x - $radius - 1],0)[1]
            + ($data[($y - 1)*$columns + $x + $radius],0)[1]
        }
      )

    return (
      for $x in 1 to $columns return $totals[$x] div ($radius * 2 + 1)
    )
  )

  return this:array($rows, $columns, $dest-data)
}

Function: vertical-blur
declare function vertical-blur($source as map(*), $radius as xs:integer) as map(*)


vertical-blur()
Quick horizontal blurring.

Params
  • source as map(*) source array
  • radius as xs:integer radius of blurring (i.e number of pixels distant)
Returns
  • map(*): array with values subject to vertical blurring
declare function this:vertical-blur (
  $source as map(*),
  $radius as xs:integer
) as map(*)
{
  this:horizontal-blur(this:transpose($source), $radius)=>this:transpose()
}

Function: blur
declare function blur($source as map(*), $radius as xs:integer) as map(*)


blur()
Quicker approximation of a Gaussian blur.

Params
  • source as map(*) source array
  • radius as xs:integer radius of blurring (i.e number of pixels distant)
Returns
  • map(*): array with values subject to Gaussian blurring
declare function this:blur(
  $source as map(*),
  $radius as xs:integer
) as map(*)
{
  this:horizontal-blur($source,$radius)=>this:vertical-blur($radius)
}

Original Source Code

xquery version "3.1";
(:~
 : 2D array of numbers: similar functionality to matrix.xqy, but with
 : sequences of doubles. Similar to the relationship between point and vector.
 : This supports more matrix operations, such as transposition and inverse.
 : This is a row major ordering: note that the get() function, like those
 : for point matrices of various kinds puts the row first e.g. get(y,x) 
 : which is, I admit a bit confusing.
 :
 : Copyright© Mary Holstege 2021-2023
 : CC-BY (https://creativecommons.org/licenses/by/4.0/)
 : @since December 2021
 : @custom:Status Stable
 :)
module namespace this="http://mathling.com/core/array"; 

import module namespace config="http://mathling.com/core/config"
       at "../core/config.xqy";
import module namespace errors="http://mathling.com/core/errors"
       at "../core/errors.xqy";
import module namespace util="http://mathling.com/core/utilities"
       at "../core/utilities.xqy";

declare namespace map="http://www.w3.org/2005/xpath-functions/map";
declare namespace math="http://www.w3.org/2005/xpath-functions/math";

declare variable $this:RESERVED := ("kind","rows","columns","data");

(:~
 : array()
 : Construct a new matrix
 : @param $rows number of rows in matrix
 : @param $columns number of columns in matrix
 : @param $data data values (must be $rows X $columns)
 : @return matrix
 : @error UTIL-BADARRAY if number of data values doesn't match size
 :)
declare function this:array(
  $rows as xs:integer,
  $columns as xs:integer,
  $data as xs:double*
) as map(*)
{
  if ($rows * $columns ne count($data))
  then errors:error("UTIL-BADARRAY", ($rows, $columns, count($data)))
  else (
    map {
      "kind": "array",
      "rows": $rows,
      "columns": $columns,
      "data": $data
    }
  )
};

(:~
 : rows()
 : Accessor for the number of rows in the array
 : 
 : @param $array the array
 : @return number of rows in array
 :)
declare function this:rows($array as map(*)) as xs:integer
{
  $array("rows")
};

(:~
 : columns()
 : Accessor for the number of rows in the array
 : 
 : @param $array the array
 : @return number of rows in array
 :)
declare function this:columns($array as map(*)) as xs:integer
{
  $array("columns")
};

(:~
 : data()
 : Accessor for the data in the array
 : 
 : @param $array the array
 : @return sequence of values (row major order)
 :)
declare function this:data($array as map(*)) as xs:double*
{
  $array("data")
};

(:~
 : data()
 : Settor for the data in the array
 : 
 : @param $array the array
 : @param $data the new data
 : @return new matrix
 : @error UTIL-BADARRAY if number of data values doesn't match size
 :)
declare function this:data($array as map(*), $data as xs:double*) as map(*)
{
  if (this:rows($array) * this:columns($array) ne count($data))
  then errors:error("UTIL-BADARRAY", (this:rows($array), this:columns($array), count($data)))
  else (
    $array=>map:put("data", $data)
  )
};

(:~
 : get()
 : Get a specific data value
 : 
 : @param $array: the array
 : @param $row: the row number
 : @param $column: the column number
 : @return data value array[row, column]
 :)
declare function this:get(
  $array as map(*),
  $row as xs:integer,
  $column as xs:integer
) as xs:double
{
  this:data($array)[($row - 1)*$array("columns") + $column]
};

(:~
 : row()
 : Return all the values from the given row
 : 
 : @param $array: the array
 : @param $row: the row number
 : @return values in row in column order
 :)
declare function this:row(
  $array as map(*),
  $row as xs:integer
) as xs:double*
{
  if ($row lt 1 or $row gt $array("rows")) then () else
  let $data := $array("data")
  let $num-columns := $array("columns")
  let $ixs :=
    for $col in 1 to $num-columns
    return ($row - 1)*$num-columns + $col
  for $ix in $ixs
  return $data[$ix]
};

(:~
 : column()
 : Return all the values from the given column
 : 
 : @param $array: the array
 : @param $column: the column number
 : @return values in column in row order
 :)
declare function this:column(
  $array as map(*),
  $column as xs:integer
) as xs:double*
{
  if ($column lt 1 or $column gt $array("columns")) then () else
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  let $ixs :=
    for $row in 1 to $num-rows
    return ($row - 1)*$num-columns + $column
  for $ix in $ixs
  return $data[$ix]
};

(:~
 : diagonal()
 : Return all the values from the diagonal
 : Raises an error if the array is not square
 : 
 : @param $array: the array
 : @return values along main diagonal
 : @error UTIL-MATRIXNOTSQUARE if matrix is not square
 :)
declare function this:diagonal($array as map(*)) as xs:double*
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    if ($num-rows != $num-columns) then (
      errors:error("UTIL-MATRIXNOTSQUARE", ($num-rows,$num-columns))
    ) else (
      let $ixs := for $i in 1 to $num-rows return ($i - 1)*$num-columns + $i
      for $ix in $ixs
      return $data[$ix]
    )
  )
};

(:~
 : trace()
 : Compute the trace of the array.
 : Raises an error if the array is not square
 : 
 : @param $array: the array
 : @return trace of the array
 : @error UTIL-MATRIXNOTSQUARE if matrix is not square
 :)
declare function this:trace($array as map(*)) as xs:double
{
  sum(this:diagonal($array))
};

(:~
 : determinant()
 : Compute the determinant of the array.
 : Raises an error if the array is not square
 : 
 : @param $array: the array
 : @return determinant
 : @error UTIL-MATRIXNOTSQUARE if matrix is not square
 :)
declare function this:determinant($array as map(*)) as xs:double
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    if ($num-rows != $num-columns) then (
      errors:error("UTIL-MATRIXNOTSQUARE", ($num-rows,$num-columns))
    ) else (
      switch($num-rows)
      case 1 return $data
      case 2 return $data[1]*$data[4] - $data[2]*$data[3]
      case 3 return (
         $data[1]*$data[5]*$data[9] + $data[2]*$data[6]*$data[7] + $data[3]*$data[4]*$data[8]
        -$data[1]*$data[6]*$data[8] - $data[2]*$data[4]*$data[9] - $data[3]*$data[5]*$data[7]
      )
      default return (
        sum(
          for $i in 1 to $num-columns return (
            math:pow(-1, 1 + $i) *
            this:determinant(
              this:array(
                $num-rows - 1,
                $num-columns - 1,
                for $r in 1 to $num-rows
                for $c in 1 to $num-columns
                where $r != 1 and $c != $i
                return $data[($r - 1)*$num-columns + $c]
              )
            )
          )
        )
      )
    )
  )
};

(:~
 : transpose()
 : Transpose the array. 
 :
 : @param $array: the array
 : @return the array, transposed
 :)
declare function this:transpose($array as map(*)) as map(*)
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    this:array(
      $num-columns,
      $num-rows,
      for $c in 1 to $num-columns
      for $r in 1 to $num-rows
      return $data[($r - 1)*$num-columns + $c]
    )
  )
};

(:~
 : cofactor-matrix()
 : Compute the cofactor matrix for the array. Raises an error if it is not square.
 :
 : @param $array: the array
 : @return matrix of cofactors
 : @error UTIL-MATRIXNOTSQUARE if matrix is not square
 : @error UTIL-CANNOT if matrix is only 1x1
 :)
declare function this:cofactor-matrix($array as map(*)) as map(*)
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    if ($num-rows != $num-columns) then (
      errors:error("UTIL-MATRIXNOTSQUARE", ($num-rows,$num-columns))
    ) else if ($num-rows = 1) then (
      errors:error("UTIL-CANNOT", ("cofactor-matrix", $data))
    ) else (
      this:array(
        $num-rows,
        $num-columns,
        for $j in 1 to $num-rows
        for $i in 1 to $num-columns
        return (
          math:pow(-1, $i + $j) *
          this:determinant(
            this:array(
              $num-rows - 1,
              $num-columns - 1,
              for $r in 1 to $num-rows
              for $c in 1 to $num-columns
              where $r != $j and $c != $i
              return $data[($r - 1)*$num-columns + $c]
            )
          )
        )
      )
    )
  )
};

(:~
 : adjoint()
 : Compute the adjoint of the array. Raises an error if it is not square.
 :
 : @param $array: the array
 : @return adjoint matrix
 : @error UTIL-MATRIXNOTSQUARE if matrix is not square
 :)
declare function this:adjoint($array as map(*)) as map(*)
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    if ($num-rows != $num-columns) then (
      errors:error("UTIL-MATRIXNOTSQUARE", ($num-rows,$num-columns))
    ) else (
      this:transpose(this:cofactor-matrix($array))
    )
  )
};

(:~
 : inverse()
 : Compute the inverse of the array. Raises an error if it is not square.
 :
 : @param $array: the array
 : @return array inverse
 : @error UTIL-MATRIXNOTSQUARE if matrix is not square
 :)
declare function this:inverse($array as map(*)) as map(*)
{
  let $data := $array("data")
  let $num-rows := $array("rows")
  let $num-columns := $array("columns")
  return (
    if ($num-rows != $num-columns) then (
      errors:error("UTIL-MATRIXNOTSQUARE", ($num-rows,$num-columns))
    ) else (
      let $det := this:determinant($array)
      let $adj := this:adjoint($array)
      return this:map($adj, function ($c as xs:double) as xs:double {$c div $det})
    )
  )
};

(:~
 : times()
 : Scalar multiplication of an array by a constant.
 :
 : @param $array: the array
 : @param $k: multiplier
 : @return k*array
 :)
declare function this:times($array as map(*), $k as xs:double) as map(*)
{
  this:map($array, function ($c as xs:double) { $c*$k })
};

(:~
 : multiply()
 : Matrix multiplication.
 :
 : @param $a1: one array
 : @param $a2: another array
 : @return a1*a2
 : @error UTIL-MATRIXNOTCOMPAT if the array sizes are incompatible
 :)
declare function this:multiply($a1 as map(*), $a2 as map(*)) as map(*)
{
  if (this:columns($a1) != this:rows($a2)) then (
    errors:error("UTIL-MATRIXNOTCOMPAT", ("multiply", this:rows($a1), this:columns($a1), this:rows($a2), this:columns($a2)))
  ) else (
    let $rows1 := this:rows($a1)
    let $columns1 := this:columns($a1)
    let $data1 := this:data($a1)
    let $rows2 := this:rows($a2)
    let $columns2 := this:columns($a2)
    let $data2 := this:data($a2)
    return (
      this:array($rows1, $columns2, 
        for $r in 1 to $rows1
        for $c in 1 to $columns2
        let $row := this:row($a1, $r)
        let $col := this:column($a2, $c)
        return sum(for $d at $i in $row return $d * $col[$i])
      )
    )
  )
};

(:~
 : add()
 : Matrix addition
 :
 : @param $a1: one array
 : @param $a2: another array
 : @return a1+a2
 : @error UTIL-MATRIXNOTCOMPAT if the array sizes are incompatible
 :)
declare function this:add($a1 as map(*), $a2 as map(*)) as map(*)
{
  if (this:columns($a1) != this:columns($a2) or this:rows($a1) != this:rows($a2)) then (
    errors:error("UTIL-MATRIXNOTCOMPAT", ("add", this:rows($a1), this:columns($a1), this:rows($a2), this:columns($a2)))
  ) else (
    let $rows1 := this:rows($a1)
    let $columns1 := this:columns($a1)
    let $data1 := this:data($a1)
    let $rows2 := this:rows($a2)
    let $columns2 := this:columns($a2)
    let $data2 := this:data($a2)
    return (
      this:array($rows1, $columns1, 
        for $d at $i in $data1 return $d + $data2[$i]
      )
    )
  )
};

(:~
 : subtract()
 : Matrix subtraction
 :
 : @param $a1: one array
 : @param $a2: another array
 : @return a1-a2
 : @error UTIL-MATRIXNOTCOMPAT if the array sizes are incompatible
 :)
declare function this:subtract($a1 as map(*), $a2 as map(*)) as map(*)
{
  if (this:columns($a1) != this:columns($a2) or this:rows($a1) != this:rows($a2)) then (
    errors:error("UTIL-MATRIXNOTCOMPAT", ("subtract", this:rows($a1), this:columns($a1), this:rows($a2), this:columns($a2)))
  ) else (
    let $rows1 := this:rows($a1)
    let $columns1 := this:columns($a1)
    let $data1 := this:data($a1)
    let $rows2 := this:rows($a2)
    let $columns2 := this:columns($a2)
    let $data2 := this:data($a2)
    return (
      this:array($rows1, $columns1, 
        for $d at $i in $data1 return $d - $data2[$i]
      )
    )
  )
};

(:~
 : map()
 : Run a function over every item in the array, return the new array, preserving properties.
 :
 : @param $array: the array
 : @param $f: function to apply to all the data in the array
 : @return new array
 :)
declare function this:map($array as map(*), $f as function(xs:double) as xs:double) as map(*)
{
  util:merge-into(
    $array,
    map {"data": for-each(this:data($array), $f)}
  )
};

(:~
 : snap()
 : Snap every value in the array to an integer. Returns an array with integer data values.
 : Preserves properties.
 :
 : @param $array: the array
 : @return new array
 :)
declare function this:snap($array as map(*)) as map(*)
{
  util:merge-into(
    $array,
    map {
      "data":
        for-each(this:data($array),
          function ($c as xs:double) as xs:integer {
            fn:round($c) cast as xs:integer
          }
        )
    }
  )
};

(:~
 : floor()
 : Snap every value in the array to the floor of that value, as an integer. 
 : Returns an array with integer data values. Preserves properties.
 :
 : @param $array: the array
 : @return new array
 :)
declare function this:floor($array as map(*)) as map(*)
{
  util:merge-into(
    $array,
    map {
      "data": 
        for-each(this:data($array), 
          function($c as xs:double) as xs:integer {
            fn:floor($c) cast as xs:integer
          }
        )
    }
  )
};

(:~
 : decimal()
 : Truncate the number of digits in every value in the array to the given number of digits, 
 : returning the new array. Preserves properties.
 :
 : @param $array: the array
 : @param $digits: number of digits to preserve
 : @return new array
 :)
declare function this:decimal($array as map(*), $digits as xs:integer) as map(*)
{
  util:merge-into(
    $array,
    map {
      "data":
        for-each(this:data($array),
          function($c as xs:double) as xs:double {
            util:decimal($c, $digits)
          }
        )
    }
  )
};

(:~
 : quote() 
 : Return a string representation of the array, suitable for debugging.
 :
 : @param $array: the array
 : @return string
 :)
declare function this:quote($array as map(*)) as xs:string
{
  string-join((
    "[",
    for $row in 1 to this:rows($array)
    return string-join(this:row($array,$row),","),
    "]"
  ),"
")
};

(:~
 : same() 
 : Equality test for two arrays. Arrays are the same if they have the same dimensions and data. 
 : Auxiliary properties do not factor into this test.
 :
 : @param $this: one array
 : @param $that: another array
 : @return true if the arrays have the same data and sizes
 :)
declare function this:same($this as map(*), $other as map(*)) as xs:boolean
{
  this:rows($this)=this:rows($other) and
  this:columns($this)=this:columns($other) and
  deep-equal(this:data($this), this:data($other))
};

(:~
 : horizontal-blur()
 : Quick horizontal blurring.
 : @param $source source array
 : @param $radius radius of blurring (i.e number of pixels distant)
 : @return array with values subject to horizontal blurring
 :)
declare function this:horizontal-blur (
  $source as map(*),
  $radius as xs:integer
) as map(*)
{
  let $rows := this:rows($source)
  let $columns := this:columns($source)
  let $data := this:data($source)

  let $dest-data := (
    for $y in 1 to this:rows($source)
    let $total1 :=
      sum(for $kx in -$radius to $radius return $data[($y - 1)*$columns + $kx])

    let $totals :=
      fold-left(2 to $columns, $total1,
        function($totals as xs:double*, $x as xs:integer) {
          $totals,
          $totals[last()]
            - ($data[($y - 1)*$columns + $x - $radius - 1],0)[1]
            + ($data[($y - 1)*$columns + $x + $radius],0)[1]
        }
      )

    return (
      for $x in 1 to $columns return $totals[$x] div ($radius * 2 + 1)
    )
  )

  return this:array($rows, $columns, $dest-data)
};

(:~
 : vertical-blur()
 : Quick horizontal blurring.
 : @param $source source array
 : @param $radius radius of blurring (i.e number of pixels distant)
 : @return array with values subject to vertical blurring
 :)
declare function this:vertical-blur (
  $source as map(*),
  $radius as xs:integer
) as map(*)
{
  this:horizontal-blur(this:transpose($source), $radius)=>this:transpose()
};

(:~
 : blur()
 : Quicker approximation of a Gaussian blur.
 : @param $source source array
 : @param $radius radius of blurring (i.e number of pixels distant)
 : @return array with values subject to Gaussian blurring
 :)
declare function this:blur(
  $source as map(*),
  $radius as xs:integer
) as map(*)
{
  this:horizontal-blur($source,$radius)=>this:vertical-blur($radius)
};