Introduction

Numeric Lua is a numerical package for the Lua programming language. It includes support for complex numbers, multidimensional matrices, random number generation, and special functions. Most of the routines are simple wrappers for stable and well known libraries from Netlib.

For instance, the special function package and the complex number package use routines from Netlib's SLATEC and DCDFLIB, the random number package uses Takuji Nishimura and Makoto Matsumoto's Mersenne Twister generator as the "engine" (uniform deviates) and Netlib's RANLIB for the remaining deviates, 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 (for a simple example, check usage below). 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.

Numeric Lua is licensed under the same license as Lua -- the MIT license -- and so can be freely used for academic and commercial purposes.

Download

Numeric Lua can be obtained from Luaforge.

The win32 binaries are provided for a quick and ready evaluation (and use); they contain a standalone version and the libraries. It should be noted that the binaries are generic, that is, they do not use any optimized version of BLAS/LAPACK.

If you got the source files, proceed to the next section for instructions.

Note: This is a very beta release, with minimal testing. Any bugs, patches, comments and suggestions are welcome.

Installation

The installation process should be simple and easy on any flavor of Unix environment (that includes MinGW and Cygwin on Windows): in a nutshell, a plain make and a further make install should do it. The process can be divided in two: basic library compilation and package compilation and installation.

After downloading the source numlua-X.Y.tar.gz, where X.Y is the current version, issue:

$ tar -zxf numlua-X.Y.tar.gz
$ cd numlua-X.Y

and edit the TARGET variable in Makefile to your platform (currently only "linux" and "win32" are available as options).

Basic library compilation

The first step is to compile the basic libraries in the lib folder. These include dcdflib, fnlib (a Netlib's SLATEC subset) and ranlib (Mersenne Twister RNG + some routines from Netlib's RANLIB). A light version of BLAS and LAPACK, containing only double and complex double routines, is also available for completeness' sake; note, however, that an optimized version of these libraries is preferred, like ATLAS or Intel's MKL.

To compile the basic libraries, run:

$ cd lib; make

from numlua-X.Y (root) dir, and wait a few minutes. :)

Some level of customization can be achieved by editing numlua-X.Y/lib/config. If you want to use an optimized version of BLAS/LAPACK, you should also comment out libblas.a and liblapack.a from LIBS in numlua-X.Y/lib/Makefile, or issue make noblas, to save some time.

Package compilation and installation

After compiling the basic libraries, change dir to src:

$ cd ../src

Edit Makefile to your needs/taste: set TARGET as in numlua-X.Y/Makefile, set your BLAS library ("blas" is set as default, in case you want to use the provided BLAS/LAPACK libraries in numlua-X.Y/lib), set installation paths, compilation flags and programs. You are now ready to compile and install:

$ make && make install

A standalone version of Numeric Lua is also available in numlua-X.Y/etc. The drill is the same: cd numlua-X.Y/etc, edit Makefile, make && make install.

Quick test

There is no testing routine yet, unfortunately. However, you can try the following commands for a quick "test drive":

$ lua
Lua 5.1.1  Copyright (C) 1994-2006 Lua.org, PUC-Rio
> -- let's generate some random numbers (fast!)
> require "numlua.rng"
> r = rng(os.time())
> for i=1,5 do print(r:unif()) end -- numbers can vary, of course :)
0.47017181466799
0.63861816644203
0.55608703091275
0.68107319541741
0.24100433581043
>
> -- now, let's add some math to spice up
> require "numlua.spfun"
> print(spfun.log1p(1e-16), math.log(1 + 1e-16)) -- keep precision
1e-16    0
> print(spfun.isinf(1/0), spfun.isnan(0/0))
true    true
> -- normal quantiles, mean = 0, var = 1
> print(spfun.qnorm(0.025, 0, 1), spfun.pnorm(spfun.qnorm(0.025)))
-1.9599639845401        0.025
>
> -- complex numbers
> require "numlua.complex"
> print(complex(1,2) + complex.i)
1 + 3i
> print(complex.exp(complex.sqrt(-2))) -- cos(sqrt(2)) + sin(sqrt(2))*i
0.15594369476537 + 0.98776594599274i
>
> -- matrices
> require "numlua.matrix"
> x = matrix(2,3,4) -- 2 x 3 x 4 matrix
> for k, i, j in x:entries() do x[k][i][j] = k * (i + j) end
> matrix.foreach(x, print) -- list x
... (contents are listed)
> -- let's try the type system
> x = matrix.zeros(1000) -- 1000 x 1000 matrix
> x:apply(function(i, j) return i + j end) -- x[i][j] = i+j
> print(x:gettype()) -- x is general
G
> -- let's profile transpose(x) * x
> local t = os.clock(); print((#x) * x); print(os.clock() - t)
matrix: 0xb586001c
1.32
> x:settype() -- guess x's type
> print(x:gettype()) -- x is symmetric
S
> local t = os.clock(); print((#x) * x); print(os.clock() - t)
1.28
> -- a bit faster... let's compute the inner product directly
> local t = os.clock(); print(x:inner()); print(os.clock() - t)
0.66

Usage

In this first example it is shown that a matrix can be used in a "standard" way, similar to a regular Lua table:

-- returns a symm posdef matrix with integer entries from Pascal's triangle
function pascal(n)
  assert(type(n) == "number", "number expected, got " .. type(n))
  assert(n > 0, "order must be positive")
  local m = matrix(n,n)
  -- first row
  for i = 1, n do m[1][i] = 1 end
  -- remaining rows
  for i = 2, n do
    m[i][1] = 1
    for j = 2, n do
      m[i][j] = m[i][j-1] + m[i-1][j]
    end
  end
  return m
end

Now, to demonstrate some functional facilities:

-- prints the entries of a matrix, right aligned
function pretty(m)
  assert(matrix.check(m), "matrix expected, got " .. type(m))
  assert(matrix.getn(m) <= 2)
  -- maximum entry length
  local ml = matrix.reduce(m, function(x, y)
    local s = string.len(y);
    if x < s then x = s end
    return x
  end, 0)
  if matrix.getn(m) == 1 then -- vector?
    matrix.foreache(m, function(e) -- for each entry
      print(string.rep(" ", 2 + ml - string.len(e)) .. e)
    end)
  else -- array
    for i = 1, m:size(1) do -- print each row
      print(matrix.reduce(m[i], function(x, y)
        return x .. string.rep(" ", 2 + ml - string.len(y)) .. y
      end, ""))
    end
  end
end

Finally, some more specialized routines:

-- create random number generator
local r = rng(os.time())
-- generate n multivariate normal deviates X_i ~ N(mean, var)
function mvnorm(mean, var, n)
  -- do some boring testing on the args ...
  -- compute square root of var using Cholesky decomposition
  local u, info = matrix.chol(var) -- #u * u = var
  assert(info == 0, "matrix is not positive definite")
  -- generate y where y:col(i) ~ N(0_m, I_m)
  local m = mean:size()
  local y = matrix(m, n):foreache(function(e) return r:norm(0, 1) end)
  -- compute #u * y + mean ~ N(mean, var)
  y = #u * y
  for i, j, e in y:entries() do
    y[i][j] = e + mean[i]
  end
  return y
end

-- solve system of linear equations A*X = B
function lsolve(A, B)
  assert(matrix.check(A), "matrix expected, got " .. type(A))
	assert(matrix.check(B), "matrix expected, got " .. type(B))
  -- try to solve by elimination first
  local x, info = matrix.linsolve(A, B)
	if info < matrix.EPS then -- A is nearly singular?
    print("WARNING: matrix is singular to machine precision: rcond = "
        .. info)
    -- try by least squares then
    x, info = matrix.lss(A, B)
  end
  return x
end

These routines can be found at lua/matrix.lua.

Reference

Unfortunately, there is no documentation available (yet!).