Introduction
What Numeric Lua isNumeric Lua is a numerical package for the Lua programming language. It includes support for complex numbers, multidimensional matrices, random number generation, fast Fourier transforms, and special functions. Most of the routines are simple wrappers for well known numerical libraries: complex numbers and part of the extended math modules come from C99; other special functions, including statistical functions, are adapted from Netlib's SLATEC and DCDFLIB; random number generation is based on Takuji Nishimura and Makoto Matsumoto's Mersenne Twister generator as the "engine" (uniform deviates) and Netlib's RANLIB for the remaining deviates; fast Fourier transforms are implemented from FFTW; and the matrix package draws most of its numeric intensive routines from Netlib's ubiquitous BLAS and LAPACK packages.
Numeric Lua tries to maintain Lua's minimalist approach by providing
bare-bone wrappers to the numerical routines. The user can use the outputs for
further computations and is then fully responsible for the results. Other Lua
features are also available, such as OO simulation through metamethods and
functional facilities. A basic API is provided in order to promote
extensibility. Also, check numlua.seeall
for a quick way to start
using Numeric Lua.
Numeric Lua is licensed under the same license as Lua -- the MIT license -- and so can be freely used for academic and commercial purposes.
numlua.buffer(opt [, arg])
Operates on Numlua buffer according to opt
:
"release"
: release all buffers with size at leastarg
;arg
defaults to 0."status"
: returns the total size and number of buffers that are busy ifarg
is true or free otherwise.
numlua.type(obj)
If obj
is an userdata with metatable mt
and
mt.__type
is not nil, returns
mt.__type
, otherwise returns the type of obj
(type(obj)
).
numlua.opmode ([mode])
Returns and optionally sets the operation mode to in-place. If
mode
is absent, returns the current operation mode; otherwise,
sets the operation mode to mode
and returns the old mode.
Note that metamethods are affected by opmode; for example, if
numlua.opmode()
is true,-a
updates
a
in-place and the statement u = u + c * v
updates
both v
and u
.
numlua.seeall
A require
call to numlua.seeall
loads most of the
functions in Numlua to the global table to make it easier to prototype and
test code or even just use Numlua is a flexible computational shell. In
particular, numlua.seeall
does the following:
- Loads the
help
module, if installed. - Sets all methods in Numlua to the global table; a prefix is
appended to a method if it exists in more than one class table or in the
global table to avoid conflicts. The order in which the methods are
"registered" and their prefixes are:
"x"
formath
(regular lua module) andmathx
,"c"
forcomplex
,"m"
formatrix
,"s"
forstat
, and"r"
forrng
. - Mathematical functions that are defined on numbers, complex numbers, and
matrices --
abs
,exp
, and so on -- are shadowed by functions that call a particular method based on its argument. For example,exp(x)
callsmath.exp(x)
ifx
is a number,complex.exp(x)
(orcexp(x)
) ifx
is complex, andmatrix.exp(x)
(ormexp(x)
) ifx
is a matrix. type
is overridden withnumlua.type
andopmode
isnumlua.opmode
.
complex
Complex numbersThe complex
module implements complex numbers according to the
C99 standard. To create a complex number a + b * i
, simply issue
complex(a, b)
or a + b * complex.i
, or even simply
a + b * i
when using numlua.seeall
. You can set the
real and imaginary parts of complex z
with z.r = x
and z.i = x
respectively. Since complex numbers are userdata,
regular assignment is a reference assignment (z = c
makes
z
a reference to c
); to actually assign values, use
z._ = c
.
Similarly to the math
library, all mathematical functions for
C99 complex numbers are available; please check your man pages for more
details about the following functions:
abs arg acos acosh asin asinh
atan atanh conj cos cosh exp
imag log logabs proj real sin
sinh sqrt tan tanh
In particular, to recover the real and imaginary parts of complex
z
, use z:real()
and z:imag()
respectively. In addition to these methods, add
,
sub
, mul
, div
, and pow
implement arithmetical methods and metamethods, that is, a + b
is
equivalent to complex.add(a, b)
and so on. The "length" of
z
, #z
, is complex.abs(z)
.
All methods accept an optional last argument -- whose default value is
given by numlua.opmode
-- that codes for in-place
applications. For instance,
> c = complex(2, -1)
> print(c:exp(), c)
3.9923240484413-6.217676312368i 2-1i
> print(c:exp(true), c) -- in-place
3.9923240484413-6.217676312368i 3.9923240484413-6.217676312368i
and while c = complex.exp(z * i) * a + b; print(c)
creates
(and discards) four temporary complex numbers,
op = numlua.opmode(true) -- in-place by default now
c = complex(0, 1); c = complex.exp(c * z) * a + b
numlua.opmode(op) -- restore previous mode
creates no temporary userdata as all operations are in-place. Note that
operand order is important even when the operation is commutative: if
numlua.opmode()
is true, a + b
is
equivalent to complex.add(a, b, true)
where the result is stored
at a
, while b + a
, or
complex.add(b, a, true)
, stores at b
.
matrix
Multidimensional matricesThe matrix
module implements multidimensional matrices. New
matrix
objects are created by specifying dimensions with
matrix.new
(uninitialized), matrix.zeros
or
matrix.ones
(filled with zeros or ones), or by directly
specifying entries in a table and using matrix
(or
matrix.fromtable
.) Matrices can be real or complex, which can be
set in the last argument of these functions. To retrieve the dimensions of a
matrix use matrix.shape
or matrix.size
; in
particular, m:size"#"
returns the number of dimensions of a
matrix m
.
To access and set entries in a matrix m
with n
dimensions you can simply use m[i1]...[in]
; for example, to set
the entry at row i
and column j
of a two-dimensional
matrix you can use m[i][j] = x
. Alternatively, you can use
matrix.get
and matrix.set
. These two functions can
also be used to produce a reference to a submatrix and to assign
whole submatrices; as an example, m:set(a)
assigns a
to m
, either element-wise if a
is a conformable
matrix or by filling m
with a
if a
is a
number or complex. You can also use matrix.copy
to copy a matrix.
More general matrix references can be created with matrix.section
or matrix.slice
.
Entries in a matrix can be set more systematically using functional
facilities: matrix.apply
sets entries based on a function that
takes indices as arguments while matrix.map
takes a function on
entries. Both functions operate in-place and in
element-order, that is, in column-major order (first indices vary
faster.) Element-order is also respected when iterating on a matrix with
matrix.entries
.
Logical methods -- matrix.find
, matrix.ifelse
,
and matrix.which
-- are helpful for filtering matrix entries.
These functions always take a condition argument, which can be: a
boolean function on each entry; or, if a number x
, the boolean
function function (e) return e == x end
; or, if a conformable
matrix x
, function (e) return e == ex end
where
ex
is the entry in x
at the same element-order index
of e
, and so, a boolean function for element-wise comparison. The
condition argument defaults to
function (e) return e ~= 0 end
.
Common mathematical operations can be carried using matrix.map
but for convenience matrix
provides the following functions:
abs acos acosh asin asinh atan
atanh cos cosh exp log sin
sinh sqrt tan tanh
All these functions take an optional second argument -- whose default value
is given by numlua.opmode
-- that specifies if the operation
should be in-place. Thus, calling m:f()
where f
is
any of the functions above, is semantically equivalent to
matrix.map(inplace and m or m:copy(), f)
.
Arithmetical methods matrix.add
, matrix.mul
,
matrix.div
, and matrix.pow
operate element-wise and
accept an optional last argument that signals in-place application to the
first argument. Matrix multiplication is implemented in three specialized
methods: matrix.mmul
for general matrix multiplication,
matrix.hemul
for Hermitian matrix multiplication, and
matrix.trmul
for triangular matrix multiplication. Arithmetical
metamethods are, not surprisingly, based on these methods and should work as
expected: __add
and __sub
are based on
matrix.add
, __mul
is based on
matrix.mul
(for scalars) and matrix.mmul
,
__mod
is matrix.ls
(left division), and
__div
is based on matrix.div
or
matrix.ls
for right division (a / b
is equivalent to
t(t(b) % t(a))
when a
and b
are
conformable matrices, where t
is matrix.t
, the
transpose operator.)
Note however, that even though metamethods make the code cleaner and more intuitive, it is often computationally more efficient to use the basic methods to avoid creation of temporary matrices. Here is a quick example:
-- updates invA to inv(A + u * v') according to the
-- Sherman-Morrison formula: inv(A + u * v') = inv(A)
-- - (inv(A) * u * v' * inv(A)) / (1 + v' * inv(A) * u)
function shermanmorrison (invA, u, v)
local x = v * invA -- v' * inv(A)
return invA:copy():mmul(invA * u, x, "n", "n", -1 / (1 + dot(x, u)))
end
The next function is equivalent,
function naiveshermanmorrison (invA, u, v)
local x = v * invA -- v' * inv(A)
return invA - invA * u * x / (1 + dot(x, u))
end
but creates two extra temporary matrices: one from
t = (invA * u) * x
and other from
t / (1 + dot(x,u))
.
While favoring in-place operations can help avoid temporary object
creation, the adoption of numlua.opmode
can make notation lighter
and less functional. It should be used with care, however, since operations
being applied to the first argument of the methods break operation
commutativity; for instance, while both matrix.add(a, b)
and
matrix.add(b, a)
return a new matrix c = a + b
,
matrix.add(a, b, true)
stores -- or better, applies -- the result
in a
and matrix.add(b, a, true)
stores the result in
b
. Another example:
function sweep (a, k)
local ck = a:col(k)
local b = ck:copy()
local ak, d = a[k], b[k]
ak:div(d, false, true) -- a[k] = a[k] / b[k], in-place
for i = 1, #a do
if i ~= k then
a[i]:add(ak, -b[i], true) -- a[i] = a[i] - b[i] * a[k], in-place
end
end
ck._ = b:div(-d, false, true) -- ck = b, b = -b / d in-place
ck[k] = 1 / d
return a
end
and now using in-place operations by default and metamethods:
numlua.opmode(true) -- in-place
function sweep (a, k)
local ck = a:col(k)
local b = ck:copy()
local ak, d = a[k], b[k]
ak = ak / d
for i = 1, #a do
if i ~= k then
a[i]:add(ak, -b[i]) -- a[i] = a[i] - ak * b[i]
end
end
ck._ = b / (-d) -- equivalent to ck:set(b / (-d))
ck[k] = 1 / d
return a
end
Finally, matrices can be serialized (in HDF5 format) using
matrix.save
and retrieved using matrix.load
.
matrix.get(m, i1 [, i2 [... [, in]]])
matrix.get(m, eoindex)
Returns a submatrix of m
. There are two ways of specifying the
submatrix:
- Using the indexes
i1
, ...,in
, wheren
ism:size"#"
, the number of dimensions inm
: returns the submatrixm[i1]...[in]
. In particular,get(m, i1, ..., in)
returns the entrym(i1, ..., in)
, andm[i]
returns thei
-th "row" ofm
. Indices that are greater than their dimension or negative are circularly remapped to the correct range: ifi
is the index for dimensiond
,i
is mapped to(i - 1) % d + 1
ifi > 0
and(i + 1) % n + n
ifi < 0
. Null indices are not allowed. - Providing an element-order vector
eoindex
; in this case, returns the respective entries inm
, that is, returnsv
such thatv[i] = m[eoindex[i]]
fori=1, ..., #eoindex
.
See also: matrix.set
matrix.set(m, i1 [, i2 [... [, in]]], value)
matrix.set(m, eoindex, value)
matrix.set(m, what, value)
Sets a submatrix of m
to value
. There are three
ways of specifying the submatrix:
- Using the indexes
i1
, ...,in
, wheren
ism:size"#"
, the number of dimensions inm
: operates on the submatrixm[i1]...[in]
. Similarly tomatrix.get
, indices can be remapped if out of range and null indices are forbidden. - Providing an element-order vector
eoindex
; in this case,m[eoindex[i]] = value[i]
, that is, the respective entries inm
are set fori=1, ..., #eoindex
. - Specifying
what
to"_"
to assign the whole matrix,"l"
or"L"
to assign only the strict lower triangle ofm
,"u"
or"U"
to assign only the strict upper triangle ofm
, and"d"
or"D"
to assign the diagonal ofm
.
Parameter value
can be a number (or complex) or a conformable
matrix according to the submatrix being indexed.
See also: matrix.get
matrix.slice(m [, first] [, last] [, step])
Returns a section of matrix m
that is a slice in the
first dimension of m
starting at first
(defaults to
1), ending at last
(defaults to #m
), and alternating
by step
(defaults to 1) positions. The returned section is a
reference to m
.
A shortcut for matrix.slice(m, first, last, step)
is
m(first, last, step)
.
See also: matrix.section
matrix.section(m, section)
Returns a section of matrix m
according to table
section
. section
is taken as a table with integer
indices from 1 to n
, the number of dimensions of m
;
the i
-th value of section
should contain another
table specifying a slice in the form {first, last,
step}
, where first
and step
default to 1
and last
defaults to m:size(i)
. The returned section
is a reference to m
.
A shortcut for matrix.section(m, section)
is
m(section)
.
Example:
-- assume a and b are real two-dimensional matrices
-- check `matrix.kronecker` for a more efficient version
function kronecker (a, b)
local ra, ca = a:shape()
local rb, cb = b:shape()
local x = matrix.new(ra * rb, ca * cb)
for i = 1, ra do
for j = 1, ca do
x{{(i-1)*rb+1, i*rb}, {(j-1)*cb+1, j*cb}}._ = b * a[i][j]
end
end
return x
end
See also: matrix.slice
matrix.add(a, b [, s [, inplace]])
Adds s * b
to a
, where s
defaults to
1. If b
is a number s
is ignored and sets
a[i] = a[i] + real(b)
if a
is real, or
a[i] = a[i] + b
if a
is complex; otherwise, sets
a[i] = a[i] + b[i] * s
if b
is a conformable matrix.
In particular, matrix.add(a, b, -1)
subtracts b
from
a
.
If inplace
is true, a
is updated
in-place.
matrix.mul(a, b [, inplace])
Performs element-wise multiplication of a
and
b
. If b
is a number, sets
a[i] = a[i] * real(b)
if a
is real, or
a[i] = a[i] * b
if a
is complex; otherwise, sets
a[i] = a[i] * b[i]
if b
is a conformable matrix.
If inplace
is true, a
is updated
in-place.
matrix.div(a, b [, right [, inplace]])
Performs element-wise division of a
and
b
. If b
is a number, sets
a[i] = a[i] / real(b)
if a
is real, or
a[i] = a[i] / b
if a
is complex; otherwise, sets
a[i] = a[i] / b[i]
if b
is a conformable matrix.
If right
is true, a
is divided
in the right, that is, a = b / a
. If inplace
is
true, a
is updated in-place.
matrix.pow(a, b [, inplace])
Performs element-wise power of a
and b
.
If b
is a number, sets a[i] = a[i] ^ real(b)
if
a
is real, or a[i] = a[i] ^ b
if a
is
complex; otherwise, sets a[i] = a[i] ^ b[i]
if b
is
a conformable matrix.
If inplace
is true, a
is updated
in-place.
matrix.new(d1, d2, ... [, iscomplex])
matrix.zeros(d1, d2, ... [, iscomplex])
matrix.ones(d1, d2, ... [, iscomplex])
Returns a matrix with dimensions d1
, d2
, and so
on. If last argument iscomplex
is true, the
returned matrix is complex. The entries are arbitrary (not
initialized) when using new
, initialized to zero when using
zeros
, and initialized to one when using ones
.
matrix.eorder(m, i1, i2, ..., in)
Returns the element-order index in m
that corresponds to the
dimensional indexes i1
, i2
, ..., in
,
where n
is the dimensionality of m
.
See also: matrix.eindex
matrix.eindex(m, eopos)
Returns the dimensional indexes in m
that correspond to
element-order index eopos
.
Example:
function list (m)
for i, e in m:entries(true) do
local t = {m:eindex(i)}
t[#t + 1] = e
print(unpack(t))
end
end
See also: matrix.eorder
,
matrix.entries
matrix.entries(m [, eorder])
Returns an iterator over the entries of m
so that the
following construction iterates over the entries of m
:
for i1, i2, ..., in, e in m:entries() do
-- body
end
Here n = m:size"#"
is the number of dimensions of
m
and e = m[i1][i2]...[in]
takes the value of each
entry in m
.
This construction is semantically equivalent to
for i = 1, m:size"*" do
local i1, i2, ..., in = m:eindex(i)
local e = m:get(i1, i2, ..., in)
-- body
end
See also: matrix.eindex
matrix.size(m [, index])
Returns the size of m
, depending on parameter
index
:
- if
index
is a number, returns the size of theindex
-th dimension ofm
- if
index == "#"
, returns the number of dimensions ofm
- if
index == "*"
, returns the total number of elements ofm
The default value for index
is 1.
See also: matrix.shape
matrix.shape(m [, index [, complex]])
Returns the dimensions of m
starting at dimension
index
(defaults to 1) and if m
is complex when
complex
is true.
See also: matrix.size
matrix.iscomplex(m)
Returns true if m
is complex and
false otherwise.
matrix.real(a)
Returns a real matrix that references the real part of
a
if a
is complex or returns a
itself
if a
is real.
See also: matrix.imag
,
matrix.complex
matrix.imag(a)
Returns a real matrix that references the imaginary part
of a
if a
is complex or returns
matrix.zeros(a:shape())
if a
is real.
See also: matrix.real
,
matrix.complex
matrix.complex(a [, b])
Returns a complex matrix m
such that
real(m) = real(a)
and imag(m) = imag(b)
. The default
value for b
is a
, so, in particular,
complex(a)
returns a copy of a
if a
is
complex. An error is raised if a
and b
are not
conformable.
See also: matrix.real
,
matrix.imag
matrix.conj(a [, inplace])
Sets each element in a
to its conjugate in-place if
inplace
is true, or returns the conjugate of
a
. If a
is not complex, does nothing.
matrix.copy(m [, x])
Returns a copy of m
and optionally sets its content to
x
, a number. matrix.copy
is semantically equivalent
to
function matrix.copy (m, x)
local c = matrix.new(matrix.shape(m, 1, true)) -- frame
return set(c, x or m)
end
See also: matrix.set
matrix.reshape(m, d1, d2, ..., dn)
Returns a reference to m
that has dimensions
d1
, d2
, ...
, dn
. The new
dimensions need to be consistent, that is,
d1 * d2 * ... * dn == m:size"*"
.
See also: matrix.size
matrix.spread(m [, dim [, count]])
Returns a replicated version of m
along dimension
d
(defaults to 1) count
number of times (defaults to
1).
matrix.find(m [, cond [, reverse]])
Returns the first element-order index in m
whose respective
entry satisfies cond
. If reverse
is
true, returns the first occurrence in m
that
does not satisfy cond
. If no entry in m
satisfies cond
, nil is returned. Note that
matrix m
is traversed only up to the first successful entry.
Examples:
function any (m, cond)
return matrix.find(m, cond) ~= nil
end
function all (m, cond)
return matrix.find(m, cond, true) == nil
end
matrix.ifelse(m, cond, x [, y])
Sets the entries in m
that satisfy cond
to
x
and, if y
is specified, the entries in
m
that do not satisfy cond
to y
.
Parameters x
and y
can be numbers or conformable
matrices to m
.
Examples:
-- note that matrix.which(m, cond, "#") is preferred
function count (m, cond)
return matrix.sum(matrix.ifelse(m:copy(), cond, 1, 0))
end
-- mask has ones for true and zeros for false
function merge (x, y, mask)
return matrix.ifelse(mask:copy(), 1, x, y)
end
local lt = function (x, y) return x < y and 1 or 0 end
function min (x, y)
return ifelse(matrix.map(x:copy(), y, lt), 1, x, y)
end
function max (x, y)
return ifelse(matrix.map(x:copy(), y, lt), 0, x, y)
end
See also: matrix.which
matrix.which(m [, cond [, what]])
Returns entries in m
that satisfy cond
depending
on what
:
what = "k"
: returns a vector of element-order indices inm
whose respective entries satisfycond
or nil otherwise;what = "v"
: returns a vector containing the entries inm
that satisfycond
or nil otherwise;what = "#"
: returns only the number of entries that safistycond
.
Parameter what
defaults to "k"
.
Examples:
-- simple alias
function count (m, cond)
return matrix.which(m, cond, "#")
end
-- mask has ones for true and zeros for false
function pack (m, mask)
return matrix.which(m, mask, "v")
end
function unpack (v, mask, m)
return matrix.set(m, matrix.which(m, mask), v)
end
See also: matrix.ifelse
matrix.apply(m, f [, eorder])
Applies a function f
to each entry in m
in-place, that is, sets m[i1,...,in] = f(i1, ..., in, m[i1,
..., in])
where i1, ..., in
are indices in
m
(m
has dimension n
). If
eorder
is true, f
takes indices in
element-order. If f
does not return a numeric value,
m
is not changed for the called index.
Examples:
local h = function (i, j) return 1 / (i + j - 1) end
function hilbert (n)
return matrix.new(n, n):apply(h)
end
local circ = function (v, n)
return function (i, j) return v[(j - i) % n + 1] end
end
function circulant (v)
local n = #v
return matrix.new(n, n):apply(circ(v, n))
end
matrix.fold(m, f [, init])
Computes the left fold of function f
over the entries
in m
in element-order, that is, returns
f(...f(f(init, m[1]), m[2])..., m[n])
See also: matrix.sum
matrix.map(m1, m2, ..., mn, f)
Sets, in element-order, the i-th entry in m1
in-place
to
m1[i] = f(m1[i], m2[i], ..., mn[i])
if the returned value is a number or complex, where m2
, ...,
mn
are all conformable to m1
.
matrix.sum(m [, alpha [, init]])
Performs a linear fold on matrix m
in element-order, that is,
computes
(...((init * alpha + m[1]) * alpha + m[2]) ... ) * alpha + m[n])
where n = m:size"*"
.
Parameters alpha
and init
default to 1 and 0
respectively; in particular, matrix.sum(m)
computes the simple
sum of the entries in m
. It is semantically equivalent to
matrix.fold(m, function(a, e) return a * alpha + e end, init)
See also: matrix.fold
matrix.min(m)
Returns the minimum element in m
and the minimizing
element-order index of m
. Complex numbers are compared based on
their real parts, and their imaginary parts are only used to break ties.
See also: matrix.max
, matrix.sort
matrix.max(m)
Returns the maximum element in m
and the maximizing
element-order index of m
. Complex numbers are compared based on
their real parts, and their imaginary parts are only used to break ties.
See also: matrix.min
, matrix.sort
matrix.sort (m [, decreasing [, returnindex]])
Sorts m
in-place. If decreasing
is
true, sorts in decreasing order. If returnindex
is true returns index
such that
m[i]
was m[index[i]]
before sorting, where
i
is in element-order; otherwise, returns m
.
As in matrix.min
and matrix.max
, complex numbers
are compared based on their real parts, and their imaginary parts are only
used to break ties.
Example:
iv = matrix.sort(v, false, true) -- increasing, return index
x = matrix.new(matrix.shape(v, 1, true))
x[iv] = v -- restore v
See also: matrix.min
, matrix.max
matrix.linspace(a, b [, n])
Returns a vector v
of length n
-- taken to be
int(|b - a| + 1)
if not specified -- such that
v[i] = (i - 1) * s
where
s
is (b - a) / (n - 1)
, and so v[1] = a
and v[n] = b
.
Examples:
v = linspace(1, n) -- v = 1:n
function arith (a, r, n)
return linspace(a, r * (n - 1) + a, n)
end
function seq (a, b, step)
local n = floor(abs((b - a) / s + 1))
return linspace(a, b, n)
end
matrix.dot(a, b [, trans])
Computes the dot product of conformable matrices a
and
b
, a^H b
as the sum of their element-wise product.
If trans
is true, the conjugate of
a
is not used (and so a^T b
is actually
computed.)
matrix.cross(a, b [, inplace])
Computes the cross product of ternary vectors a
and
b
in-place in a
if inplace
is
true.
matrix.col(m, i)
Returns a reference to the i
-th column of the
two-dimensional matrix m
.
matrix.transpose (m [, hermitian])
Returns the transpose of matrix m
. If hermitian
is true, returns the conjugate transpose of
m
.
matrix.diag(m [, shift])
If m
is a vector, returns a square matrix such that the
shift
diagonal of the matrix is m
; otherwise, if
m
is two-dimensional, returns a reference to the
shift
diagonal of m
. If shift
is zero
(the default value), the main diagonal is referenced; if
shift > 0
, the shift
diagonal above the
main diagonal is referenced, and if shift < 0
, the
shift
diagonal below the main diagonal is referenced.
Example:
function trace (m)
local r, c = m:shape()
assert(m:size"#" == 2 and r == c, "square matrix expected")
return m:diag():sum()
end
matrix.concat(m1, ..., mn [, colcat])
Concatenates conformable matrices m1
, ...
,
mn
row-wise. If colcat
is
true, concatenates column-wise.
matrix.c(v1, ..., vn)
Concatenates vectors or numbers v1
, ...
,
vn
into a new vector.
matrix.swap (a, b)
Swaps the contents of conformable vectors a
and
b
.
matrix.pivot(m, pvt [, colpivot [, inplace]])
Given pivot vector pvt
, returns a new matrix p
with the same rows of m
but with rows i
and
pvt[i]
swapped. Note that, in particular, p = P * m
,
where P = eye(#m):pivot(pvt)
, that is, P
is the
permutation matrix for pvt
.
If colpivot
is true, columns are swapped
instead of rows, that is, returns p = m * P
where
P = eye(#m):pivot(pvt, true)
.
If inplace
is true, rows/columns are swapped
in-place.
Examples:
function perm (p) -- column permutation
local n = #p
local s = matrix.linspace(1, n):pivot(p)
local P = matrix.zeros(n, n)
for i = 1, n do P[i][s[i]] = 1 end
return P
end
-- see `matrix.lu` and `matrix.qr` for details:
require "numlua.seeall"
l, u, p = lu(a)
print(pretty(perm(p) * a)) -- = a:pivot(p)) = l * u
q, r, p = qr(a,true)
print(pretty(q * r * t(perm(p)))) -- = (q * r):pivot(p,true) = a
matrix.trmul (x, A [, uplo [, invert [, trans [, side]]]])
For conformable two-dimensional matrices A
and x
,
A
square, returns the result of tri(A) * x
in-place in x
, where tri(A)
takes the lower
triangle section of A
if uplo == "l"
or
uplo == "L"
, or the upper triangle section of A
if
uplo == "u"
or uplo == "U"
. The default value for
uplo
is "L"
.
If invert
is true, multiplies by the inverse
of A
.
If trans == "t"
or trans == "T"
, returns
tri(A)^T * x
; if trans == "c"
or
trans == "C"
, returns tri(A)^H * x
(conjugate
transpose); if trans == "n"
or trans == "N"
, does
nothing (default behavior).
If side == "l"
or side == "L"
, A
is
multiplied at the left of x
(default); if
side == "r"
or side == "R"
, A
is
multiplied at the right, that is, returns x * tri(A)
.
For examples, see matrix.lu
and matrix.qr
.
matrix.hemul (C, A, [, inner, [, what [, alpha]]])
Given a Hermitian matrix C
, performs the Hermitian update
C = C + alpha * A^H * A
if inner
is
true, or C = C + alpha * A * A^H
otherwise. The
default value for alpha
is 1.
Parameter what
specifies how C
is to be updated;
if what == "l"
or what == "L"
, only the lower
triangle of C
is updated; if what == "u"
or
what == "U"
, only the upper triangle of C
is
updated; if what == "f"
or what == "F"
, the lower
triangle of C
is updated by the multiplication and reflected to
the strict upper triangle. The default value for what
is "F".
Example:
-- generates Wishart deviates
function wishart (S, n)
local m = #S
local L = assert(chol(S, "L")) -- Cholesky factor: S = L * L^T
local c = zeros(m) -- cache
return function ()
local w = zeros(m, m)
for i = 1, n - 1 do
c = trmul(rnorm(0, 1, c), L) -- c ~ N(0, S)
w = hemul(w, c, false, "L") -- w = w + c * c^T
end
return hemul(w, trmul(rnorm(0, 1, c), L)) -- n-th run: full hemul
end
end
matrix.mmul (C, A, B [, transA, transB, alpha])
For conformable matrices A
, B
, and
C
, updates C
with the multiplication
C = C + alpha * transA(A) * transB(B)
. Parameters
transA
and transB
specify possible operations for
A
and B
:
"n"
or"N"
does nothing (default);"t"
or"T"
takes the transpose;"c"
or"C"
takes the conjugate transpose.
If B
is a vector, transB
is ignored. Moreover,
A
can only be a vector if B
is also a vector, in
which case the outer product A * B^H
is computed regardless of
transA
and transB
.
matrix.norm(m [, what])
Computes the norm of m
according to parameter what
:
- if
what
is a numberp
, computes the Lp norm what == "m"
orwhat == "M"
: computesmax(abs(m))
what == "o"
orwhat == "O"
: computes the one norm, that is, the maximum column sum ifm
is a two-dimensional matrixwhat == "i"
orwhat == "I"
: computes the infinity (sup) norm, that is, the maximum row sum ifm
is a two-dimensional matrixwhat == "f"
orwhat == "F"
: computes the Frobenius norm
The default value for what
is "F"
. If
what == "m"
, what == "M"
,
what == "i"
, or what == "I"
, also returns the index
the maximizes the norm.
matrix.chol(A [, what [, inplace]])
Computes the Cholesky factorization of symmetric positive-definite matrix
A
. Parameter what
defines the computation:
- If
what == "l"
orwhat == "L"
, the default value, references only the lower triangle ofA
and returns a lower triangular matrixL
such thatA = L * L^H
; - If
what == "u"
orwhat == "U"
, references only the upper triangle ofA
and returns an upper triangular matrixU
such thatA = U^H * U
.
If an internal error occurs, returns nil and an error
message; if A
is not positive-definite, returns
false and a warning message. If inplace
is
true, overwrites A
with its Cholesky factor.
Examples:
function matrix.isposdef (m)
local c, msg = chol(m)
if c == nil then error(msg) end
return not c == false
end
-- generates multivariable normal deviates (also check rng.rmvnorm)
function mvnorm (mu, S)
local L = assert(chol(S, "L")) -- Cholesky factor: S = L * L^T
return function (dest)
local dest = dest or zeros(#mu)
local s = rnorm(0, 1, dest) -- s ~ N(0, I_n)
s = trmul(s, L) -- s = L * s, s ~ N(0, S)
s = add(s, mu, true) -- s = s + mu, s ~ N(mu, S)
return s
end
end
See also: matrix.hemul
,
matrix.trmul
matrix.lu(A [, inplace])
Computes the LU decomposition of A
: returns an unit lower
triangular matrix L
, an upper triangular matrix U
and a permutation vector pvt
such that
A:pivot(pvt) = L * U
, or equivalently,
A = P^T * L * U
where P = eye(#A):pivot(pvt)
. If
inplace
is true, both L
and
U
are returned in-place at A
,
L
on the strict lower triangle and U
on the
upper triangle.
Examples:
-- compute determinant of m
function det (m)
local c = assert(lu(copy(m), true)) -- in-place on copy
return prod(diag(c))
end
-- solve linear system A * x = b
function linsolve (A, b)
local l, u, p = lu(A)
local x = b:pivot(p) -- x = eye(#A):pivot(p) * b
x:trmul(l, 'l', true) -- x = inv(l) * x
x:trmul(u, 'u', true) -- x = inv(u) * x
return x
end
See also: matrix.pivot
matrix.rcond(A [, what [, inplace]])
Estimates the reciprocal of the condition number of matrix A
in 1-norm, depending on parameter what
:
what == "d"
orwhat == "D"
, assumes thatA
is a diagonal matrix;what == "l"
orwhat == "L"
, assumes thatA
is a lower triangular matrix, ifwhat == "u"
orwhat == "U"
, assumesA
is an upper triangular matrix;what == "p"
orwhat == "P"
, assumes thatA
is a Hermitian positive-definite matrix;what == "g"
orwhat == "G"
, the default value, assumes thatA
is a general matrix.
If inplace
is true and what
codes for A
as being Hermitian positive-definite or general,
A
is overwritten with either a Cholesky or LU decomposition
respectively.
See also: matrix.svd
, matrix.inv
matrix.inv(A [, what [, inplace [, norcond]]])
Returns the inverse of square matrix A
and an estimate of the
reciprocal of its condition number, according to parameter
what
:
what == "d"
orwhat == "D"
, assumes thatA
is a diagonal matrix;what == "l"
orwhat == "L"
, assumes thatA
is a lower triangular matrix, ifwhat == "u"
orwhat == "U"
, assumesA
is an upper triangular matrix;what == "p"
orwhat == "P"
, assumes thatA
is a Hermitian positive-definite matrix;what == "g"
orwhat == "G"
, the default value, assumes thatA
is a general matrix.
If inplace
is true and what
codes for A
as being Hermitian positive-definite or general,
A
is overwritten with either a Cholesky or LU decomposition
respectively.
If norcond
is true, the reciprocal of the
condition number is not estimated.
Example:
-- Square-root of a matrix A by Denman-Beavers algorithm
function dbsqrtm (A, tol, maxiters)
local opmode = numlua.opmode(true) -- set in-place operations
local iY, iZ = zeros(A:shape()), zeros(A:shape()) -- buffers
local s, n = inf, 1 -- norm(Y), #iterations
local Y, Z = copy(A), eye(#A) -- initialize
while true do
iY._, iZ._ = Y, Z
Y, Z = (Y + inv(iZ)) / 2, (Z + inv(iY)) / 2 -- in-place
-- check termination:
local f = norm(Y)
if abs(f - s) <= tol or n > maxiters then break end
s, n = f, n + 1
end
numlua.opmode(opmode) -- restore previous mode
return Y, Z, s
end
See also: matrix.rcond
matrix.svd(A [, what])
Computes the singular value decomposition of
m
-by-n
matrix A
. The SVD is specified
as A = U * S * V^H
, where U
and V
are
unitary matrices whose first min(m, n)
columns represent the left
and right singular vectors of A
, respectively, and S
is a diagonal matrix whose min(m, n)
diagonal entries contain the
singular values of A
. U
, V
, and
S
are computed according to parameter what
:
what == "n"
orwhat == "N"
, returns only the singular values ofA
;what == "l"
orwhat == "L"
, the firstmin(m, n)
columns ofU
are stored in-place atA
;what == "r"
orwhat == "R"
, the firstmin(m, n)
rows ofV^H
are stored in-place atA
;what == "a"
orwhat == "A"
, the default value, returnsU
,S
, andV^H
.
The singular values are returned in descending order. Note that
V^H
is returned, not V
.
If an internal error occurs, returns nil and a diagnostic message; if the algorithm fails to converge, returns false and a warning message.
Examples:
function cond (m)
local s = assert(svd(m, "n")) -- just singular values
return s[1] / s[#m] -- condition number
end
local lt = function (x) return function(e) return e < x end end
function matrix.rank (m, tol) -- effective rank according to tol
local s = assert(svd(m, "n")) -- just singular values
local tol = tol or 0
if tol <= 0 then -- set default tolerance?
tol = max(m:shape()) * eps * s[1]
end
local r = s:find(lt(tol))
if r == nil then -- all > tol?
return #s -- full rank: min(m:shape())
end
return r - 1 -- last position > tol
end
For other examples, see matrix.null
, matrix.orth
,
and matrix.pinv
in matrix.lua
.
See also: matrix.rcond
matrix.qr(A [, permute [, inplace]])
Computes the QR factorization of matrix A
, that is, returns an
unitary matrix Q
and an upper trapezoidal matrix R
such that A = Q * R
.
If permute
is true, computes the QR
factorization with column pivoting by also returning a permutation vector
p
such that A = (Q * R):pivot(p, true)
, that is,
P = eye(A:size(2)):pivot(p, true)
satisfies
A * P^T = Q * R
.
In inplace
is true, A
is
overwritten with R
.
Example:
-- solve linear system A * x = b
function linsolve (A, b)
local q, r = qr(A)
local x = zeros(#b):mmul(A, b, 't') -- x = t(A) * b
x:trmul(r, 'u', true, 't') -- x = t(inv(r)) * x
x:trmul(r, 'u', true) -- x = inv(r) * x
return x
end
matrix.eig(A [, what [, hermitian]])
Computes the eigenvalues of A
, and optionally the eigenvectors
according to parameter what
:
what == "n"
orwhat == "N"
: do not compute the eigenvectors ofA
, just return the eigenvalues ofA
(in a vector);what == "l"
orwhat == "L"
: returns the eigenvalues (in a vector) and the left eigenvectors ofA
(in a matrix);what == "r"
orwhat == "R"
: the default option, returns the eigenvalues (in a vector) and the right eigenvectors ofA
(in a matrix);what == "a"
orwhat == "A"
: returns the eigenvalues and both left and right eigenvectors ofA
.
Note that the right eigenvectors of A
are the left eigenvalues
of A^H
and vice-versa. If hermitian
is
true, A
is assumed to be Hermitian and only its
upper triangle is accessed.
When computing the eigenvalues -- and optionally eigenvectors -- of
A
, the matrix is balanced to improve accuracy; see
matrix.balance
for a description of the operations. There is
currently no way of preventing balancing.
Examples:
-- condition number
function cond (m)
local x = zeros(#m, #m):hemul(m) -- x = m * t(m)
local e = eig(x, "n", true) -- sqrt(e) are singular values of m
return math.sqrt(math.abs(max(e) / min(e)))
end
-- if F[n] is the n-th Fibonacci number,
-- {{F[n-1], F[n]}, {F[n], F[n+1]}} = {{0, 1}, {1, 1}} ^ n
local s, V = eig(matrix{{0, 1}, {1, 1}}, "r", true)
local v = V[1] * V[2] -- element-wise multiplication
local w = zeros(2) -- workspace
function fib (n)
w._ = s; w:pow(n, true) -- w = s ^ n, in-place
return dot(v, w) -- V * diag(s ^ n) * t(V)
end
See also: matrix.balance
,
matrix.svd
matrix.balance (A [, what])
Balances a matrix A
as a way to reduce its 1-norm and thus
improve the accuracy of the computed eigenvalues and eigenvectors. The
operations involve, first, permuting A
by a similarity
transformation to isolate eigenvalues in the first 1 to l
and
last h
to #A
elements of the diagonal, and then
scaling by a diagonal similarity transformation to rows and columns
l + 1
to h - 1
to make them as close in norm as
possible. Parameter what
controls which operations are performed
to balance A
:
what == "n"
orwhat == "N"
: nothing is done (no balancing);what == "p"
orwhat == "P"
: permute only;what == "s"
orwhat == "S"
: scale only;what == "b"
orwhat == "B"
: both permute and scale, the default option
See also: matrix.eig
matrix.ls (A, b [, svd [, tol [, inplace]]])
Computes the least-squares solution to the linear system
A * x = b
, that is, finds x
that minimizes the
2-norm of b - A * x
, where A
is a
m
-by-n
matrix. Parameter b
can be a
(m
-by-nrhs
) matrix containing many right-hand side
vectors as columns.
If svd
is true, the solution is found using a
singular value decomposition; otherwise, a QR factorization is applied. The
effective rank of A
is relative to tol
, which
defaults to max(m, n) * mathx.eps
(see Example in
matrix.svd
for details).
If inplace
is true, the solution is returned
in b
if the system is not underdetermined, that is, if
A
has at least as many rows as columns. In this case,
b:slice(1,n)
contains the solutions x
and the sum of
squares for the i
-th column of b:slice(n+1, m)
is
the residual sum of squares for the solution of the i
-th
right-hand side.
The solution and the effective rank are returned; in addition, if
svd
is true, the singular values of
A
are also returned (as a vector).
Example:
-- basic LS linear model fitting
function matrix.lm (a, b, svd)
local m, n = a:shape()
assert(m >= n, "system is underdetermined")
assert(b:size"#" == 1, "single RHS expected")
local x, rank = ls(a, b, svd)
-- report summary statistics
local coef = x:slice(1, n)
local rss = (b - a * coef):norm() ^ 2
local rss0 = (b - b:sum() / m):norm() ^ 2
local df = m - rank
local F = df / (rank - 1) * (rss0 / rss - 1)
local pvalue = 1 - stat.pf(F, rank - 1, df)
return {coef = coef, rss = rss, df = df, F = F, pvalue = pvalue}
end
See also: matrix.svd
, matrix.qr
matrix.fft(m [, inverse [, wisdomonly [, inplace]]])
Computes the fast Fourier transform of complex matrix
m
. If inverse
is true, computes the
normalized inverse FFT. wisdomonly
is a FFTW-specific
parameter: if true, a plan to compute the transform is only
created if wisdom is available for the given problem, that is, a
FFTW_WISDOM_ONLY
flag is passed when building the plan. The
default flag is FFTW_ESTIMATE
. If inplace
is
true, m
is updated in-place.
See also: matrix.fct
, fft.plan
matrix.fct(m [, inverse [, wisdomonly [, inplace]]])
Computes the discrete cosine transform of real matrix
m
. If inverse
is true, computes the
normalized inverse DCT. wisdomonly
is a FFTW-specific
parameter: if true, a plan to compute the transform is only
created if wisdom is available for the given problem, that is, a
FFTW_WISDOM_ONLY
flag is passed when building the plan. The
default flag is FFTW_ESTIMATE
. If inplace
is
true, m
is updated in-place.
See also: matrix.fft
, fft.plan
matrix.save(m, filename)
Saves matrix m
in HDF5 file filename
.
matrix.load(filename)
Loads a matrix from HDF5 file filename
and returns the loaded
matrix.
fft
Fast Fourier transforms
The fft
module implements fast Fourier transforms. The simple
interface to FFTs comprehends just matrix.fft
and
matrix.fct
which operate directly on a matrix. If, however, you
intend to execute the same transform many times it might be more efficient to
use a plan, so check fft.plan
. If more control on how
plans are executed is needed, there is also an interface to the
wisdom of FFTW in fft.wisdom
.
Example:
> a = matrix({{1, 2, 3}, {0, -1, -4}}, true)
> b = matrix(a:shape(1, true))
> p = fft.plan(b, false, fft.flag["measure"]) -- direct transform
> ip = fft.plan(b, true, fft.flag["measure"]) -- inverse transform
> b._ = a
> p(); print(b:pretty()) -- execute direct
1+0i 1-1.7320508075689i 1+1.7320508075689i
11+0i -4+3.4641016151378i -4-3.4641016151378i
> ip(); print(b:pretty()) -- execute inverse
1+0i 2+0i 3+0i
0+0i -1+0i -4+0i
fft.plan(m [[, inverse], flag])
Returns a FFTW plan that uses matrix m
as a buffer to apply
discrete transforms in-place in m
. If m
is
complex, a fast Fourier transform is executed, otherwise, if m
is
real, a discrete cosine transform is performed. If inverse
is
true, normalized inverse transforms are
computed.
Parameter flag
specifies a planner flag. Possible values are
available in table fft.flag
. The default flag is
fft.flag["measure"]
, FFTW_MEASURE
. Note that only
"estimate"
(FFTW_ESTIMATE
) and
"wisdomonly"
(FFTW_WISDOM_ONLY
) do not
overwrite buffer m
; for all other flags, m
has to be
initialized after the plan is created.
If plan = fft.plan(m, ...)
, #plan
returns
m
and plan()
executes the plan.
See also: matrix.fft
,
matrix.fct
fft.wisdom([wisdom])
Provides an interface for FFTW wisdoms. If wisdom
is
nil, returns the currently accumulated wisdom; if
wisdom
is true, forgets accumulated wisdom; if
wisdom
is a string, imports wisdom from it.
stat
Statistics and probability methods
The stat
module implements statistical routines: probability
mass and density functions, cumulative distribution functions (cdf) and
inverse cdf (quantile) functions for a number of distributions; and
factors.
stat.factor(t)
Constructs a factor based on table t
, a
sequence. A factor is a more compact representation of t
where each distinct entry -- a level -- is assigned an integer key,
the key level. If f
is a factor built from
t
, f[i]
returns t[i]
, f(i)
returns the key level of t[i]
, and f()
returns the
levels of f
as a table; in particular, f[i]
is
equivalent to f()[f(i)]
. The expression #f
returns a
vector with level counts, that is, (#f)[l]
returns the number of
entries in t
with key level l
and thus values
f()[l]
; note that (#f):sum()
is #t
.
Example:
> f = stat.factor{"a", "c", "a", "a", "b", "c"}
> for l, level in ipairs(f()) do print(l, level, (#f)[l]) end
1 a 3
2 c 2
3 b 1
> for i = 1, (#f):sum() do print(i, f[i], f(i)) end
1 a 1
2 c 2
3 a 1
4 a 1
5 b 3
6 c 2
See also: factor.fold
,
factor.partition
, factor.design
factor.fold(f, m, func [, init])
Folds the entries in matrix m
according to the levels in
factor f
and fold function func
, that is, for each
key level l
, returns a vector v
such that v[l]
= func(m[i1], func(m[i2], func(..., func(m[in], init)...)))
where
i1, ..., in
satisfy f(i1) = ... = f(in) = l
. The
default value for init
is zero.
Example:
> f = stat.factor{"a", "c", "a", "a", "b", "c"}
> add2 = function(x, y) return x + y end
> x = matrix.linspace(1, (#f):sum())
> for l, v in f:fold(x, add2):entries() do print(l, v) end
1 8
2 8
3 5
See also: stat.factor
,
factor.partition
factor.partition(f, m)
Partitions the matrix m
according to the levels in
f
by returning a table t
such that t[l]
is a vector containing the entries i1, ..., in
in m
such that f(i1) = ... = f(in) = l
.
Example:
> f = stat.factor{"a", "c", "a", "a", "b", "c"}
> x = matrix.linspace(1, (#f):sum())
> vprint = function(v) return table.concat(v:totable(), ", ") end
> for l, v in ipairs(f:partition(x)) do print(l, vprint(v)) end
1 1, 3, 4
2 2, 6
3 5
See also: stat.factor
,
factor.fold
factor.design(f, [ref])
Retuns design matrix for factor f
. If ref
is zero
(default value) then no contrasts are assumed; otherwise, return a "sum to
zero" contrast matrix with ref
as reference level.
Example:
> f = stat.factor{"a", "c", "a", "a", "b", "c"}
> print(f:design():pretty())
1 0 0
0 1 0
1 0 0
1 0 0
0 0 1
0 1 0
> print(f:design(2):pretty())
1 1 0
-1 1 -1
1 1 0
1 1 0
0 1 1
-1 1 -1
stat.dbeta(x, a, b)
stat.pbeta(x, a, b)
stat.qbeta(p, a, b)
Computes the density (dbeta
), cumulative distribution function
(pbeta
), and inverse cdf (qbeta
) for the
beta distribution with parameters a >= 0
and
b >= 0
. The density is given by
x^(a - 1) * (1 - x)^(b - 1) / B(a, b)
for
0 <= x <= 1
, where B(a, b)
is the beta
function.
See also: mathx.beta
stat.dbinom(s, xn, pr)
stat.pbinom(s, xn, pr)
stat.qbinom(p, xn, pr)
Computes the probability mass function (dbinom
), cumulative
distribution function (pbinom
), and inverse cdf
(qbinom
) for the binomial distribution with size
xn >= 0
and probability of success on each trial
0 <= pr <= 1
. The density is given by
choose(xn, s) * pr^s * (1 - pr)^(xn - s)
for
s = 0, ..., xn
, where choose(xn, s)
is the binomial
coefficient.
stat.dchisq(x, df [, pnonc])
stat.pchisq(x, df [, pnonc])
stat.qchisq(p, df [, pnonc])
Computes the density (dchisq
), cumulative distribution
function (pchisq
), and inverse cdf (qchisq
) for the
chi-squared distribution with df > 0
degrees of
freedom and non-centrality parameter pnonc
. If
pnonc ~= 0
, df
can be zero.
stat.dexp(x [, r])
stat.pexp(x [, r])
stat.qexp(p [, r])
Computes the density (dexp
), cumulative distribution function
(pexp
), and inverse cdf (qexp
) for the
exponential distribution with rate r
. The default value
for r
is one. The exponential distribution has density
r * exp(-r * x)
.
stat.df(x, dfn, dfd)
stat.pf(x, dfn, dfd [, pnonc])
stat.qf(p, dfn, dfd [, pnonc])
Computes the density (df
), cumulative distribution function
(pf
), and inverse cdf (qf
) for the F
distribution with dfn >= 0
and dfd >= 0
degrees of freedom in the numerator and denominator, and non-centrality
parameter pnonc
.
stat.dgamma(x, shape [, rate])
stat.pgamma(x, shape [, rate])
stat.qgamma(p, shape [, rate])
Computes the density (dgamma
), cumulative distribution
function (pgamma
), and inverse cdf (qgamma
) for the
gamma distribution with parameters shape >= 0
and
rate >= 0
. The default value for rate
is one. A
gamma distribution with shape s
and rate r
has
density r^s / G(s) * x^(s - 1) * exp(-r * x)
where G
is the gamma function, for x >= 0
.
See also: mathx.gamma
stat.dhyper(x, r, b, n)
stat.phyper(x, r, b, n)
Computes the probability mass function (dhyper
) and cumulative
distribution function (phyper
) for the hypergeometric
distribution with parameters r, b >= 0
and
0 <= n <= r + b
. The density is given by
choose(r, x) * choose(b, n - x) / choose(r + b, n)
for
x = 0, ..., n
.
stat.dnbinom(s, xn, pr)
stat.pnbinom(s, xn, pr)
stat.qnbinom(p, xn, pr)
Computes the probability mass function (dnbinom
), cumulative
distribution function (pnbinom
), and inverse cdf
(qnbinom
) for the negative binomial distribution with
size xn >= 0
and probability of success on each trial
0 <= pr <= 1
. The density is given by
choose(xn + s - 1, s) * pr^xn * (1 - pr)^s
for
s = 0, ..., xn
, where choose(xn + s - 1, s)
is the
binomial coefficient.
stat.dnorm(x [, mean [, sd]])
stat.pnorm(x [, mean [, sd]])
stat.qnorm(p [, mean [, sd]])
Computes the density (dnorm
), cumulative distribution function
(pnorm
), and inverse cdf (qnorm
) for the
normal distribution with parameters mean
and standard
deviation sd >= 0
. The default values for mean
and sd
are zero and one, respectively.
stat.dpois(s, l)
stat.ppois(s, l)
stat.qpois(p, l)
Computes the probability mass function (dpois
), cumulative
distribution function (ppois
), and inverse cdf
(qpois
) for the Poisson distribution with mean
l >= 0
. The density is given by
l^s * exp(-l) / s!
for s = 0, 1, ...
.
stat.dt(x, df)
stat.pt(x, df)
stat.qt(p, df)
Computes the density (dt
), cumulative distribution function
(pt
), and inverse cdf (qt
) for the Student
t distribution with df >= 0
degrees of freedom.
rng
Random number generation and deviatesThe rng
module provides random number generators
based on the Mersenne-Twister generator and netlib's ranlib
.
Objects can be created with rng.new
, copied with
rng.copy
, and have their seed set with rng.seed
.
The following methods sample deviates from a number of probability distributions or elements in a matrix object:
rbeta rchisq rexp rf rgamma rnorm
rmvnorm runif runifx rdirichlet rbinom rnbinom
rpois runifint sample lsample
Each one of these methods are closed on a rng
object and each
object r
has a (uservalue) table containing its methods, that is,
methods for which r
is an upvalue. To retrieve the method table
for r
, use #r
. For convenience, a "class"
rng
object c
is always created and has its methods
stored in the class table so that rng.f
is
(#c).f
.
All random deviate methods accept an optional last parameter v
that must be a real matrix that serves as destination: if
v
is provided, v:size()
deviates are generated and
stored in-place at v
.
Examples:
> n = 5
> for i = 1, n do print(i, rng.rnorm()) end
1 -0.83608755044204
2 0.30498104946162
3 -1.227283848856
4 -0.91148894025705
5 0.28245551259444
> r = rng.new()
> v = (#r).rnorm(0, 1, matrix.new(n))
> for i = 1, n do print(i, v[i]) end
1 -0.83608755044204
2 0.30498104946162
3 -1.227283848856
4 -0.91148894025705
5 0.28245551259444
rng.new([seed])
Returns a new rng object r
with seed seed
. To
obtain a table with methods that use r
, call #r
.
See also: rng
, rng.seed
rng.seed([seed])
Sets the seed of rng
to seed
; seed
can be a number or a real vector.
rng.rbeta(a, b [, dest])
Returns random deviates from the beta distribution with shape parameters
a
(a > 0
) and b
(b > 0
).
See also: stat.dbeta
rng.rchisq(df [,xnonc [, dest]])
Returns random deviates from the chi-square distribution with degrees of
freedom df
(df > 0
) and optional non-centrality
parameter xnonc
(xnonc >= 0
, defaults to
zero).
See also: stat.dchisq
rng.rexp([av [, dest]])
Returns random deviates from the exponential distribution with mean
av
(defaults to one).
See also: stat.dexp
rng.rf(dfn, dfd [, xnonc [, dest]])
Returns random deviates from the F distribution with degrees of freedom
dfn
and dfd
(dfn > 0
and
dfd > 0
) and optional non-centrality parameter
xnonc
(xnonc >= 0
, defaults to zero).
See also: stat.df
rng.rgamma(a [, s [, dest]])
Returns random deviates from the gamma distribution with shape
a
and optional rate s
.
See also: stat.dgamma
rng.rnorm([m [, sd [, dest]]])
Returns random deviates from the normal distribution with mean
m
(defaults to zero) and standard deviation sd
(sd > 0
, defaults to one).
See also: stat.dnorm
, rng.rmvnorm
rng.rmvnorm(m, S [, dest])
Returns random deviates from the multivariate normal distribution with mean
m
, a real vector of size n
, and "standard deviation"
S
-- either a square matrix of order n
, taken as the
lower triangular Cholesky decomposition of the variance-covariance matrix, or
a positive real vector of size n
representing the square root of
the diagonal of the variance-covariance matrix.
Returns S * u + m
, where u
is sampled from
N(0, I_n)
.
See also: matrix.chol
, rng.rnorm
rng.runif([low [, high [, dest]]])
Returns random deviates from the uniform distribution between
low
(defaults to zero) and high
(defaults to
one).
See also: stat.dunif
rng.runifx([low [, high [, dest]]])
Returns random deviates from the uniform distribution between
low
(defaults to zero) and high
(defaults to one)
with 53-bit resolution.
See also: stat.dunif
rng.rdirichlet(alpha [, dest])
Returns random deviates from the Dirichlet distribution with parameter
alpha
. alpha
should be a real vector with positive
entries.
rng.rbinom(n, p [, dest])
Returns random deviates from the binomial distribution with size
n
(n >= 0
) and probability p
(0 <= p <= 1
).
See also: stat.dbinom
rng.rnbinom(n, p [, dest])
Returns random deviates from the negative binomial distribution with size
n
(n >= 0
) and probability p
(0 <= p <= 1
).
See also: stat.dnbinom
rng.rpois(l [, dest])
Returns random deviates from the Poisson distribution with mean
l >= 0
.
See also: stat.dpois
rng.runifint([low [, high [, dest]]])
Returns random deviates from the discrete uniform distribution between
low
and high
.
rng.sample(m [, normalized])
Returns a sample from 1 to m:size"*"
with weights specified by
the entries in the matrix m
. If normalized
is
true, the entries of m
are considered to be
normalized, that is, it is assumed that m:sum() == 1
.
rng.lsample(m [, normalized])
Returns a sample from 1 to m:size"*"
with log-weights
specified by the entries in the matrix m
. If
normalized
is true, the entries of
m
are considered to be normalized, that is, it is assumed that
m:exp():sum() == 1
.
mathx
Extended math libraryThis module contains C99 math functions, Bessel functions (from Netlib's AMOS library), and other utilitary functions. C99 math functions are listed below:
isfinite
,isinf
,isnan
,isnormal
, andsignbit
take one number argument and return a boolean;acosh
,asinh
,atanh
,cbrt
,erf
,erfc
,exp2
,expm1
,lgamma
,log1p
,log2
,logb
,nearbyint
,rint
,round
,tgamma
, andtrunc
take one number argument and return a number;copysign
,fdim
,fmax
,fmin
,hypot
,nextafter
, andremainder
take two arguments and return a number;fpclassify
takes a number and returns a string;fma
takes three numbers and returns a number;scalbn
takes a number and an "integer" and returns a number
Please refer to the man pages for details about each of these functions.
mathx.airya (z, deriv, scaling)
Returns the complex Airy function Ai(z) or its derivative if
deriv
is true, and a completion
status
. If scaling
is true, returns
Ai(z) * exp(zeta) where zeta = 2 / 3 * z * sqrt(z) to remove
the exponential decay in -pi / 3 < arg(z) < pi / 3 and the
exponential growth in pi / 3 < abs(arg(z)) < pi.
If status
is true, the return is normal, that
is, the computation completed without errors, underflows, or overflows;
otherwise, status
is an error message. If no computation is
performed, nil is returned.
See also: mathx.airyb
mathx.airyb (z, deriv, scaling)
Returns the complex Airy function Bi(z) or its derivative if
deriv
is true, and a completion
status
. If scaling
is true, returns
Bi(z) * exp(-abs(real(zeta))) where zeta = 2 / 3 * z *
sqrt(z), to remove the exponential behavior in both the left and right
half planes.
If status
is true, the return is normal, that
is, the computation completed without errors, underflows, or overflows;
otherwise, status
is an error message. If no computation is
performed, nil is returned.
See also: mathx.airya
mathx.besselh(nu, z, secondkind, scaling [, n])
Returns a n
member sequence of complex Hankel (Bessel)
functions CY(j) = H(m, nu + j - 1, z) where m = 2 if
secondkind
is true or m = 1 otherwise,
with nonnegative orders nu + j - 1 for j = 1, ..., n and
complex z ~= 0 in the cut plane -pi < arg(z) < pi. The
size of the sequence n
defaults to 1. A completion
status
is also returned.
If scaling
is true, returns CY(j) *
exp(-mm * z * I), where mm = 3 - 2 * m and I is the
imaginary unit. The scaling removes the exponential behavior in both the upper
and lower half planes.
If status
is true, the return is normal, that
is, the computation completed without errors, underflows, or overflows;
otherwise, status
is an error message. If no computation is
performed, nil is returned.
See also: mathx.besseli
,
mathx.besselj
, mathx.besselk
,
mathx.bessely
.
mathx.besseli(nu, z, scaling [, n])
Returns a n
member sequence of complex modified Bessel
functions of first kind CY(j) = I(nu + j - 1, z) with nonnegative
orders nu + j - 1 for j = 1, ..., n and complex
z
in the cut plane -pi < arg(z) < pi.
The size of the sequence n
defaults to 1. A completion
status
is also returned.
If scaling
is true, returns CY(j) *
exp(-abs(real(z))) to remove the exponential growth in both the left
and right half planes.
If status
is true, the return is normal, that
is, the computation completed without errors, underflows, or overflows;
otherwise, status
is an error message. If no computation is
performed, nil is returned.
See also: mathx.besselh
,
mathx.besselj
, mathx.besselk
,
mathx.bessely
.
mathx.besselj(nu, z, scaling [, n])
Returns a n
member sequence of complex Bessel functions of
first kind CY(j) = J(nu + j - 1, z) with nonnegative orders nu +
j - 1 for j = 1, ..., n and complex z
in the cut
plane -pi < arg(z) < pi. The size of the sequence
n
defaults to 1. A completion status
is also
returned.
If scaling
is true, returns CY(j) *
exp(-abs(imag(z))) to remove the exponential growth in both the lower
and upper half planes.
If status
is true, the return is normal, that
is, the computation completed without errors, underflows, or overflows;
otherwise, status
is an error message. If no computation is
performed, nil is returned.
See also: mathx.besselh
,
mathx.besseli
, mathx.besselk
,
mathx.bessely
.
mathx.besselk(nu, z, scaling [, n])
Returns a n
member sequence of complex modified Bessel
functions of second kind CY(j) = K(nu + j - 1, z) with nonnegative
orders nu + j - 1 for j = 1, ..., n and complex z ~=
0
in the cut plane -pi < arg(z) < pi. The size of the
sequence n
defaults to 1. A completion status
is
also returned.
If scaling
is true, returns CY(j) *
exp(z) to remove the exponential behavior in both the left and right
half planes.
If status
is true, the return is normal, that
is, the computation completed without errors, underflows, or overflows;
otherwise, status
is an error message. If no computation is
performed, nil is returned.
See also: mathx.besselh
,
mathx.besseli
, mathx.besselj
,
mathx.bessely
.
mathx.bessely(nu, z, scaling [, n])
Returns a n
member sequence of complex Bessel functions of
second kind CY(j) = Y(nu + j - 1, z) with nonnegative orders nu +
j - 1 for j = 1, ..., n and complex z
in the cut
plane -pi < arg(z) < pi. The size of the sequence
n
defaults to 1. A completion status
is also
returned.
If scaling
is true, returns CY(j) *
exp(-abs(imag(z))) to remove the exponential growth in both the lower
and upper half planes.
If status
is true, the return is normal, that
is, the computation completed without errors, underflows, or overflows;
otherwise, status
is an error message. If no computation is
performed, nil is returned.
See also: mathx.besselh
,
mathx.besseli
, mathx.besselj
,
mathx.besselk
.
mathx.feq(x1, x2 [, epsilon])
Approximate equal up to precision epsilon
, based on Knuth's
suggestion in: Knuth, D. E. (1998), "The Art of Computer Programming. Vol. 2:
Seminumerical Algorithms", Section 4.2.2.
Note that if x1 ~== x2
when feq(x1, x2)
is
true, then similarly
x1 ~< x2
whennot feq(x1, x2) and x1 < x2
x1 ~> x2
whennot feq(x1, x2) and x1 > x2
x1 ~<= x2
whenfeq(x1, x2) or x1 < x2
x1 ~>= x2
whenfeq(x1, x2) or x1 > x2
mathx.log1pe(x)
Computes log(1 + exp(x))
by avoiding calls to log
and exp
when abs(x) < eps
.
mathx.lse(x1, x2)
Computes the "soft max" of x1
and x2
, the log sum
of exps log(exp(x1) + exp(x2))
.
mathx.beta(x, y)
mathx.lbeta(x, y)
Computes the beta function B
with arguments x
and
y
and its logarithm (beta
and lbeta
respectively.) The beta function is defined by
B(x, y) = integral_0^1 t^(x - 1) (1 - t)^(y - 1) dt
, or,
equivalently, by B(x, y) = G(x) * G(y) / G(x + y)
where
G
is the gamma function.
See also: mathx.gamma
,
mathx.lgamma
mathx.digamma(x)
Computes the digamma function, the logarithmic derivative of the gamma
function, at x
.
See also: mathx.gamma
mathx.choose(n, k)
mathx.lchoose(n, k)
Compute binomial coefficients and their logarithms (choose
and
lchoose
respectively`.)
License
Copyright (c) 2011 Luis Carvalho
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.


Get the most out of Stable Diffusion 2.1: Tips and Tricks
The latest version of Stable Diffusion, with its increased support for larger 768px images, is a significant improvement over previous models that only supported 512px images. While it may require some adjustment as it was trained on a more limited dataset and prompts in a different way, this new model has the potential to create stunning AI-generated artwork. In this series, we will provide some helpful tips and techniques for utilizing the capabilities of the Stable Diffusion 2.1 model to its full potential.
Requirements
To get started with Stable Diffusion, you’ll need to install a few tools. As each of these tools comes with its own installation instructions, be sure to refer to the manuals for specific installation instructions for each tool. Here, we’ll just highlight settings specifically required in this workflow. As Stable Diffusion web UI needs to be installed though Git, make sure you have that available. For Windows, I prefer to install Git SCM over the official client.
- Stable Diffusion web UI from AUTOMATIC1111. The UI we’ll use to generate images, train models, …
- Dreambooth Extension for Stable Diffusion web UI. Can be installed through the UI (Extensions tab)
- Stable Diffusion 2.1 768 model. Get
v2-1_768-ema-pruned.ckpt
from huggingface.co and put this file in the correct folder - Lama Cleaner AI based in-painting tool, useful for removing parts of an image (optional)
- ChaiNNer workflow editor to create image editing workflows for e.g. up-scaling (optional)
- GIMP Image editing tool, free and open source (optional)
- You’ll need up-to-date drivers for your NVIDIA GPU along with cuDNN
To get the most out of the techniques in this series, you’ll want to make sure you have a top-notch GPU. If you don’t have one with at least 16GB VRAM at your disposal, you can make use of Runpod.io to rent a machine in the cloud with a powerful GPU for less than half a dollar per hour. Specs of my machine, used throughout this series are:
- NVIDIA 4080 RTX (16 Gb VRAM)
- RAM: 32 Gb
- CPU: Ryzen 7 3700X
- 100 Gb of free space
- Windows 10
Though for this first post a more modest GPU will suffice, any NVIDIA GPU with 6Gb VRAM should work. Only when training new styles or adding persons/objects, which will be covered later, more VRAM is required.
Configuring Stable Diffusion web UI
After installing Stable Diffusion web UI, you’ll need to make a few adjustments to the file webui-user.bat
.
The program needs to be started with two additional arguments --xformers
and --no-half
, to do this add these
to the line which starts with set COMMANDLINE_ARGS=
. It should be set COMMANDLINE_ARGS=--xformers --no-half
.
--xformers
is an optimization that shaves off some VRAM usage. --no-half
is required by SD2.1.
Now save the file and run it to start Stable Diffusion (note that the first time it will install some additional
dependencies).
Once started go to Stable Diffusion web UI with your browser (default address is http://localhost:7860/) and go to the
Settings
tab. Find the option Use cross attention optimizations while training
and make sure it is
checked. This enables xformers to be used while training an embedding or hypernetwork as well, this will be covered
in a future post.
Creating art with AI
Personally I don’t think of myself as an artist, when working with AI generated art you assume the role of an art director. You provide instructions, aka. a prompt, to the neural network and hopefully tickle the right neurons to get an output you like, using a few rounds of revisions to get closer and closer to an image that meets your expectations. The workflow that works for me starts with prompt engineering, generating a large set of images, selecting the best one, generating some variations, picking the best one, use in-painting/img2img/compositing to get to the details right and finally up-scaling the image.
Step 1: Prompt Engineering
Writing a prompt is giving the neural network an indication what you need.
an steam engine leaving the trainstation, fall, sunset, painting, fine-art, detailed
For SD2.1 negative prompts are quite important, below is a set of terms I start with. In case I want a more painterly style I’ll add photo, photorealistic here (as that is what I don’t want). In case certain objects start popping up that are not wanted this would be the place to write that down and push the AI away from generating images containing those objects.
disfigured, kitsch, ugly, oversaturated, grain, low-res, Deformed, blurry, bad anatomy, disfigured, poorly drawn face,
mutation, mutated, extra limb, ugly, poorly drawn hands, missing limb, blurry, floating limbs, disconnected limbs,
malformed hands, blur, out of focus, long neck, long body, ugly, disgusting, poorly drawn, childish, mutilated,
mangled, old, surreal
After a few tries I got this image, I like the style and the locomotive reminds me of a model train I had as a kid. So let’s continue with this one.
Stable Diffusion also give you the prompts and settings used:
an steam engine leaving the trainstation, fall, sunset, painting, fine-art, detailed
Negative prompt: photo, photorealistic, disfigured, kitsch, ugly, oversaturated, grain, low-res, Deformed, blurry, bad anatomy, disfigured, poorly drawn face, mutation, mutated, extra limb, ugly, poorly drawn hands, missing limb, blurry, floating limbs, disconnected limbs, malformed hands, blur, out of focus, long neck, long body, ugly, disgusting, poorly drawn, childish, mutilated, mangled, old, surreal
Steps: 30, Sampler: Euler a, CFG scale: 7, Seed:<PHONE_NUMBER>, Size: 768x768, Model hash: 4bdfc29c, Batch size: 4, Batch pos: 3
It will take some back-and-forth changing settings to find a good prompt with setting that yield a good initial image. You can also go brute-force here and generate a dozens, or even hundreds, of images and pick the best one.
Step 2: Generating Variations
Next, we’ll generate some variations of this image to see if we can tease out a better version. To do this copy the
seed to the corresponding box and check Extra
this will show the variation options. Set variation strength to
0.05 - 0.1 (relatively low as we don’t want to change the image too much) and start generating more images. Don’t change
the prompt or the settings.
After generating a few batches the following image appeared. The building here is much closer to a train station. While it is far from perfect it is a step in the good direction. So let’s go with this one!
an steam engine leaving the trainstation, fall, sunset, painting, fine-art, detailed
Negative prompt: photo, photorealistic, disfigured, kitsch, ugly, oversaturated, grain, low-res, Deformed, blurry, bad anatomy, disfigured, poorly drawn face, mutation, mutated, extra limb, ugly, poorly drawn hands, missing limb, blurry, floating limbs, disconnected limbs, malformed hands, blur, out of focus, long neck, long body, ugly, disgusting, poorly drawn, childish, mutilated, mangled, old, surreal
Steps: 30, Sampler: Euler a, CFG scale: 7, Seed:<PHONE_NUMBER>, Size: 768x768, Model hash: 4bdfc29c, Batch size: 8, Batch pos: 1, Variation seed:<PHONE_NUMBER>, Variation seed strength: 0.05
Note that this doesn’t work with all samplers, with Euler a it seems to work quite well, however, some have
issues reproducing an image from the same seed when not generating an entire batch of images. (A solution here is to
build several batches of one image) Efficient samplers, which converge
after few steps, also generate wildly different images even with low variation (< 0.1). Hence, your milage might vary!
Another way to get different versions of the initial image is to vary the CFG Scale
and the number of
Sampling Steps
. Using the X/Y Plot script (which can be enabled at the very bottom of the text2img settings),
all combinations can be generated. Alternatively, you can skip through to step 4 and use img2img with high denoising
setting to generate variations.
Step 3: In-painting to correct mistakes
The rails in the above image are obviously incorrect, there is an extra rail we’ll need to instruct the AI to remove.
This is exactly what in-painting is designed to do, you paint over a specific part, provide a new prompt to cover that
up. Playing with the settings and being patient is important here. I masked the rail I didn’t want, and entered the
prompt ballast, pebbles, painting, fine-art, detailed
as this is what I wanted to replace the rail with.
I tried with and without Inpaint at full resolution
and ultimately landed on 0.9
for the Denoising
strength
giving the AI a lot of freedom to change the image. For adding completely new content or removing something
setting Masked content
to Latent noise
or Latent nothing
might help generating something not based
on the underlying image. In case you wish to remove something you can also move over to Lama Cleaner which is
excellent to quickly get rid of some unwanted items quickly.
Now our image looks like this:
Removing the rail in the bottom corner using in-painting turned out to be much more cumbersome. Here I loaded the
image in GIMP and used the clone tool to copy over some grass. Next, using in-painting with a low Denoising
strength
of 0.2 with a prompt of grass and dirt, some more variation in texture was added in again in this part of
the image.
Step 4: Img2Img to fine-tune the details
Now we move our image over to img2img to fine tune the details. Enter the original prompt again, set the Denoising
strength
to a low value e.g. 0.1 to 0.4 to prevent the AI from changing too much and generate a few more images.
You’ll see that different iterations will have different aspects changed, some versions have better wheels, others
have better rails, more detail in the building, …
The tick here is to pick each image that has the best something (you can quickly drag and drop images to your Desktop). Load the images as different layers in GIMP, add a layer mask and paint in the parts of the image you wish to keep. Since there is very little difference between the images this should be very easy, and you don’t have to be super precise while painting.
We’ll export the image from GIMP as we are getting close to the final version!
Step 3 and 4 can be repeated if needed.
Step 5: Up-scaling
SD2.1 allows us to start with an image that is 768x768 or larger, and hence has twice as many pixels as the 512x512 models, so up-scaling is considerably easier. Probably, sending the image over to Extras and selecting two methods will work nicely to bump the image size up 2x - 4x.
Stable Diffusion web UI also supports SD up-scaling in the img2img tab. This will chop up the image in smaller parts,
enlarge each, fill in details using img2img and combine them again. Make sure the Denoising strength
is set to
a small value. Result here vary as img2img will also alter parts, even at low denoising values, and might undo our
earlier work in-painting and compositing.
Another option is to use ChaiNNer and create an up-scale workflow, this supports additional models for up-scaling which can be combined with other filters to improve the clarity or sharpness of the image. CupScale which also supports various neural up-scalers can be used here as well.
Step 6: Final cleanup
There are still a few spot where the image has items I don’t like, for instance the ornaments in front, the antenna in the back, … As I don’t want to replace them, simple remove them Lama Cleaner is easier and faster than SD in-painting (using the lama model).
Conclusion
AI-generated art has the power to transform the way we create and view art. Our final image, reminiscent of the vintage box-art of model kits, is a testament to this. Without the help of AI tools, generating an image like this would have been impossible for me. However, while simply writing a prompt can lead to great results, it usually takes a bit more effort to get the image you desire. That’s why it’s important to keep experimenting withs setting and refining your images. Don’t be afraid to generate and re-generate your image multiple times until you’re satisfied with the outcome.
Liked this post ? You can buy me a coffee


COVID-19 Tracker
Select State
--Choose State--
District
Total Recovered
Total Confirmed Cases
Total Active


SUM Function
The SUM function returns the sum of the input column or expression values. The SUM function works with numeric values and ignores NULL values.
Syntax
SUM ( [ DISTINCT | ALL ] expression )
Arguments
- expression
-
The target column or expression that the function operates on.
- DISTINCT | ALL
-
With the argument DISTINCT, the function eliminates all duplicate values from the specified expression before calculating the sum. With the argument ALL, the function retains all duplicate values from the expression for calculating the sum. ALL is the default.
Data Types
The argument types supported by the SUM function are SMALLINT, INTEGER, BIGINT, NUMERIC, DECIMAL, REAL, and DOUBLE PRECISION.
The return types supported by the SUM function are
-
BIGINT for BIGINT, SMALLINT, and INTEGER arguments
-
NUMERIC for NUMERIC arguments
-
DOUBLE PRECISION for floating point arguments
The default precision for a SUM function result with a 64-bit NUMERIC or DECIMAL argument is 19. The default precision for a result with a 128-bit NUMERIC or DECIMAL argument is 38. The scale of the result is the same as the scale of the argument. For example, a SUM of a DEC(5,2) column returns a DEC(19,2) data type.
Examples
Find the sum of all commissions paid from the SALES table:
select sum(commission) from sales; sum ------------- 16614814.65 (1 row)
Find the number of seats in all venues in the state of Florida:
select sum(venueseats) from venue where venuestate = 'FL'; sum -------- 250411 (1 row)
Find the number of seats sold in May:
select sum(qtysold) from sales, date where sales.dateid = date.dateid and date.month = 'MAY'; sum ------- 32291 (1 row)


Jurgen Klopp states why Alex Oxlade-Chamberlain has not been starting more games
Alex Oxlade-Chamberlain has only started twice in the Premier League this season for Liverpool as he continues to recover from his knee injury
by Alex SmithJurgen Klopp has explained that he has kept Alex Oxlade-Chamberlain out of his Liverpool starting XI because he is still finding his rhythm.
Midfielder Oxlade-Chamberlain has started two Premier League games this season and has also kicked off against Liverpool and MK Dons in the Carabao Cup, along with the 4-1 victory over Genk in the Champions League.
The England international only managed 19 minutes of action last season following his horrendous knee injury and despite putting in some eye-catching performances this season, with three goals in his last two appearances, is yet to return to full steam - by his own admission.
Oxlade-Chamberlain said earlier this week: "Quality-wise, I am not too happy with my performance to be honest.
"I wasn't great on the ball and not near my own personal level of where I want to be. I put in a good shift and run a lot.
"I watch the boys of a weekend and I know what it is about here. You have got to be at a really high physical level to go and put in a performance like that.
And Klopp has admitted that that is music to his ears and has deemed Oxlade-Chamberlain a "long-term project".
"He has scored, nice, both at Genk, a real stunner on Wednesday, a lot of outstanding moments, but other moments he is not involved.
"We always have to try to reach perfection. He still needs rhythm but that is no problem. He is a long-term project.
"It is good a player thinks he is not at his best. He could think he is because he scored some goals. Quite some space for improvement. That is good.
"Best shape ever before injury? He was exceptional."


{% include meta.html %}
{{ site.name }}
{{ site.description }}
Parents+Riders
Race+Travel
News
Coaches
Sponsors
Forms+Links
{{ content }}
{% include analytics.html %}


First Name*
Last Name
UserName*
Email*
Password*
Re-enter Password*
Enter DOB:
You are:
Qualifications
Nationality


Marathon by the Sea
Saint John NB/CAN
Race Winners
(23 runnings)
23
Sep
2007
13
Aug
2006
07
Aug
2005
08
Aug
2004
10 Aug 2003
11 Aug 2002
19
Aug
2001
20
Aug
2000
1999
09 Aug 1998
17 Aug
1997
18 Aug 1996
1995
Notes: 1995 was the inaugural year.
13
Aug 2017
2:37:09
Matthew McNeil (NS/CAN)
3:34:53
Charlotte Mawer-2
14
Aug 2016
2:48:44
Ryan O'Shea-3
3:40:20
Charlotte Mawer (NB/CAN)
09 Aug
2015
2:42:25
Ryan O'Shea-2
3:25:26
Andree Germain (PQ/CAN)
10 Aug 2014
2:51:38
Ryan O'Shea (NB/CAN)
3:02:53
Shelley Doucet (NB/CAN)
11
Aug 2013
2:52:10
Tim Harte (PA/USA)
3:27:12
Elita Rahn (NS/CAN)
12
Aug 2012
2:45:16
Chris Brake-3
3:28:24
Lisa Leskien (ON/CAN)
14
Aug 2011
2:37:19
Chris Brake-2
3:27:16
Cara Allaway (ON/CAN)
08 Aug
2010
2:48:08
Alex Coffin-5
3:30:49
Margaret Johnson-2
27 Sep
2009
2:40:01
Chris Brake (NB/CAN)
3:25:45
Brenda Guitard (NB/CAN)
28 Sep
2008
2:51:00
Leo McCosham (PE/CAN)
3:28:00
Margaret Johnson (NB/CAN)
2:50:42
Alex Coffin-4
3:05:10
Heather Goodfellow-3
2:41:55
Dave Eaton (NS/CAN)
3:16:13
Tammy Mercier (NS/CAN)
2:47:05
David MacLennan-3
3:32:17
Allyson MacDonald (NB/CAN)
2:48:39
Alex Coffin-3
3:13:12
Heather Goodfellow-2
2:40:20
David MacLennan-2
3:12:18
Heather Goodfellow (NB/CAN)
2:35:56
David MacLennan (NS/CAN)
3:26:45
Shelley Gauthier-2
2:55:14
Alex Coffin-2
3:12:31
Rebecca Richards (NS/CAN)
2:41:56.1
Victor Hickey (ON/CAN)
3:13:56.8
Shelley Gauthier (NB/CAN)
---
n/a
---
n/a
2:35:02
Alex Coffin (NB/CAN)
3:08:57
Mary Bartlett (NB/CAN)
---
n/a
2:58:11
Donna Lockett (NY/USA)
---
n/a
---
n/a
---
n/a
---
n/a
Most Wins (1995/date)
Men
5Alex Coffin (NB/CAN) 1998,2001,2004,2007,2010
Women
3Heather Goodfellow (NB/CAN) 2003-04,2007
Prize Money (ARRS qualified
runners only)
No information on prize money has been
received for this race.
Race Time Bias
(sec/km)
Insufficient data to calculate the race time bias.
Competitive Level
(1990-date)
Insufficient data to calculate competitive levels.
Quality Performances
Men No men under 2:20 Women No women under 2:50
Last Updated on 14 Aug 2017
Contributor: Maurice Wilson (5)