External package for storing the exported symbols.

Internal package for implementing numcl.

NUMCL defines several symbols which have the same name as the corresponding CL symbols. We call them conflicting symbols. To avoid the confusion in the code base, we use 3 packages: NUMCL.EXPORTED, NUMCL.IMPL, NUMCL.

This package replaces COMMON-LISP package by shadowing-import symbols from NUMCL.EXPORTED on top of COMMON-LISP package.

Returns a string that represents the version of NUMCL currently installed.

Returns true when ARRAY satisfies the NUMCL assumption, that is, an array displaced to a non-displaced 1D array.

Type specifier for the arrays satifying the NUMCL assumption, that is, an array displaced to a non-displaced 1D array.

Returns a simple array of the equivalent contents.

Reshape the array while sharing the backing 1D array. -1 implies that the axis size is deduced from the other axes. At most one axis is allowed to be -1. T implies that the axis size is preserved. It can be used as many times, but only at the right/leftmost axes.

Example of reshaping (3 8 5):

valid:

(6 -1 10)     = (6 2 10)
(t 2 2 2 t)   = (3 2 2 2 5)
(3 t t)       = (3 8 5)
(2 -1 2 2 t)  = (2 3 2 2 5)

invalid:

(2 t 2 2 t)

axes: an int or a list of ints

An extended aref that accepts ranges as lists, similar to numpy's array access. For a 3D array x,

• range

x[1:5,2,3]   = (aref x '(1 5) 2 3)
x[2,1:5,3]   = (aref x 2 '(1 5) 3)
x[2,1:2:5,3] = (aref x 2 '(1 2 5) 3)
x[2,1:,3]    = (aref x 2 '(1 t) 3)
x[2,:1,3]    = (aref x 2 '(t 1) 3)
x[2,:,3]     = (aref x 2 '(t t) 3)
x[2,:,3]     = (aref x 2    t   3)

• insufficient axis

(aref x '(1 5)) == (aref x '(1 5) t t)
(aref x 2 '(1 5)) == (aref x 2 '(1 5) t)

• newaxis

(aref x '(1 2 5) nil 2 3)

• ellipsis

(aref x '- 2) = (aref x t t 2) = x[...,2]
(aref x 2 '-) = (aref x 2 t t) = x[2,...]
(aref x 2 '- 3) = (aref x 2 t 3) = x[2,...,3]
(aref x 2 3 '-) = (aref x 2 3 t) = x[2,3,...]

Arange's argument signature is irregular, following the API of numpy. The interpretation of its arguments depends on the number of arguments.

(arange stop &key type) (arange start stop &key type) (arange start stop step &key type)

Don't worry, we provide a compiler-macro to avoid the runtime dispatch.

Note: ASARRAY is SLOW as it recurses into the substructures.

Copy CONTENTS to a new array. When CONTENTS is a multidimentional array, its elements are copied to a new array that guarantees the NUMCL assumption. When CONTENTS is a nested sequence, it is traversed up to the depth that guarantees the sane shape for an array. When elements are copied, it is coerced to TYPE. When TYPE is not given, it is replaced with the float-contagion type deduced from the elements of CONTENTS. It may return a 0-dimensional array with CONTENTS being the only element.

For example:

;; a vector of two lists.
(asarray '((1) (1 2)))               -> #((1) (1 2))
;; a 2D array of 4 lists.
(asarray '(((1) (1 2)) ((3) (3 4)))) -> #2A(((1) (1 2)) ((3) (3 4)))

(asarray '((1 2) (3 4)))    -> #2A((1 2) (3 4))
(asarray #(#(1 2) #(3 4)))  -> #2A((1 2) (3 4))
(asarray #((1 2) (3 4)))    -> #2A((1 2) (3 4))

However, this behavior may not be ideal because the resulting shape could be affected by the lengths of the strings.

(asarray #(#(1 2) #(3 4)))   -> #2A((1 2) (3 4))
(asarray #(#(1 2) #(3 4 5))) -> #(#(1 2) #(3 4 5))

(asarray #("aa" "aa"))   -> #2A((#a #a) (#a #a))
(asarray #("aa" "aaa"))  -> #("aa" "aaa")

As a remedy to this problem, we allow TYPE to be a specifier for vector subtypes. Providing such a type specifier will keep the leaf objects (e.g. strings) from split into individual elements. We don't allow it to be a multidimentional array [at the moment.]

(asarray #(#(1 2) #(3 4))   :type '(array fixnum (*))) -> #(#(1 2) #(3 4))
(asarray #(#(1 2) #(3 4 5)) :type '(array fixnum (*))) -> #(#(1 2) #(3 4 5))

(asarray #("aa" "aa")  :type 'string)    -> #("aa" "aa")
(asarray #("aa" "aaa") :type 'string)    -> #("aa" "aaa")

(asarray '((1 2) (3 4))   :type '(array fixnum (* *)))  -> error

Note: Broadcasting with EINSUM is experimental.

Performs Einstein's summation. The SUBSCRIPT specification is significantly extended from that of Numpy and can be seens as a full-brown DSL for array operations.

SUBSCRIPTS is a sequence of the form (<ISPEC>+ [-> <TRANSFORM>*] [-> <OSPEC>*] [-> <IOPTION>*] [-> <OOPTION>*]). The remaining arguments ARGS contain the input arrays and optionally the output arrays.

Inputs are interpreted in the following rule, depending on the number of arrows.

 (<ISPEC>+)
 (<ISPEC>+ -> <OSPEC>*)
 (<ISPEC>+ -> <TRANSFORM>* -> <OSPEC>*)
 (<ISPEC>+ -> <TRANSFORM>* -> <OSPEC>* -> <IOPTION>*)
 (<ISPEC>+ -> <TRANSFORM>* -> <OSPEC>* -> <IOPTION>* -> <OOPTION>*)

SPEC

The first set of SPECs specifies the input subscripts, and the second set of SPECs specifies the output subscripts. Unlike Numpy, there can be multiple output subscripts: It can performs multiple operations in the same loop, then return multiple values. The symbol -> can be a string and can belong to any package because it is compared by STRING=.

Each SPEC is an string designator, such as a symbol IJK or a string "IJK", where each character is considered as an index. Broadcasting can be specified by a hyphen. It signals a type-error when it contains any non-alpha char except hyphens (-). Note that a symbol NIL is interpreted as an empty list rather than N, I and L. Once broadcasting is used in one input, it must be used in all inputs/outputs.

Alternatively, each SPEC can be a list that contains a list of symbols. For example, ((i j) (j k) -> (i k)) and (ij jk -> ik) are equivalent.

When -> and the output SPECs are omitted, a single output is assumed and its spec is a union of the input specs. For example, (ij jk) is equivalent to (ij jk -> ijk). Note that (ij jk) and (ij jk ->) have the different meanings: The latter sums up all elements.

TRANSFORM

TRANSFORM is a list of element-wise operations. The number of TRANSFORM should correspond to the number of outputs. In each TRANSFORM, the elements in the input arrays can be referenced by $N, where N is a 1-indexed number. Similarly the output array can be referred to by @N.

• Example: (ij ik -> (+ @1 (* $1 $2)) -> ik) is equivalent to (ij ik -> ik) (a GEMM).
• Example: (i -> (sin $1) -> i) is equivalent to mapping a sin function to an array.

By default, TRANSFORM is (+ @1 (* $1 ... $N)) for N inputs, which is equivalent to Einstein's summation.

OPTION

Each OPTION is a list of iteration specifiers. Iteration specifier is a form (index &key (start 0) (end -1) (step 1)). START and END are the forms specifying the interval designator of the loop range, which supersedes the default full-width loop. STEP is a form specifying the increment of the index. INDEX is a symbol which should match one of the subscripts used in one of the SPEC.

For example, you can specify

 (einsum '(ij jk -> (+ @1 (* $1 $2)) -> ik -> ((j :step 2)) ((j :step 3)))
         (ones '(5 10))
         (ones '(15 5)))

to suggest that the index j for the first array should have a stride of 2, and the index j for the second array should have a stride of 3.

This is fairly complex, but we hope to provide the maximum flexibility.

ARGS

The shape of each input array should unify against the corresponding input spec. For example, with a spec IJI, the input array should be of rank 3 as well as the 1st and the 3rd dimension of the input array should be the same. Note that this is affected by OPTIONS --- when specified, it uses the number of iteration instead of the size of the dimension.

The shape of each output array is determined by the corresponding output spec. For example, if SUBSCRIPTS is (ij jk -> ik) and the input arrays are NxM and MxL matrices, the output has a shape NxL.

If the output arrays are provided, their shapes and types are also checked against the corresponding output spec. The types should match the result of the numerical operations on the elements of the input arrays.

The outputs are calculated in the following rule.

• The output array types are calculated based on the TRANSFORM, and the shapes are calcurated based on the SPEC and the input arrays.

• The output arrays are allocated and initialized by zeros.

• Einsum nests one loop for each index in the input specs. For example, (ij jk -> ik) results in a triple loop.

• In the innermost loop, each array element is bound to $1..$N / @1..@N.

• For each @i, i-th transform is evaluated and assigned to @i.

• If the same index appears multiple times in a single spec, they share the same value in each iteration. For example, (ii -> i) returns the diagonal elements of the matrix.

When TRANSFORMs are missing, it follows naturally from the default TRANSFORM values that

• When an index used in the input spec is missing in the output spec, the axis is aggregated over the iteration by summation.
• If the same index appears across the different input specs, the element values from the multiple input arrays are aggregated by multiplication. For example, (ij jk -> ik) will perform (setf (aref a2 i k) (* (aref a0 i j) (aref a1 j k))) when a0, a1 are the input arrays and a2 is the output array.

For example, (einsum '(ij jk) a b) is equivalent to:

 (dotimes (i <max> <output>)
   (dotimes (j <max>)
     (dotimes (k <max>)
       (setf (aref <output> i j k) (* (aref a i j) (aref b j k))))))

Performance

If SUBSCRIPTS is a constant, the compiler macro builds an iterator function and make them inlined. Otherwise, a new function is made in each call to einsum, resulting in a large bottleneck. (It could be memoized in the future.)

The nesting order of the loops are automatically decided based on the specs. The order affects the memory access pattern and therefore the performance due to the access locality. For example, when writing a GEMM which accesses three matrices by (setf (aref output i j) (* (aref a i k) (aref b k j))), it is well known that ikj-loop is the fastest among other loops, e.g. ijk-loop. EINSUM reorders the indices so that it maximizes the cache locality.

Equivalent of the same function in numpy. Note the default type difference.

• empty : does not explicitly fill the array. In an unsafe compiler setting, junk value may appear.
• full : fill the array with a certain value.
• zeros, ones : fill the array with zeros / ones. type affects the actual value being filled.
• X-like : similar to above functions, but takes another array and returns the array of the same shape.

Reverses the axes of an array.

Matrix product of two arrays.

Dot product of two vectors. For complex values, the first value is conjugated.

Inner product of two vectors.

Compute the outer product of two vectors.

Compute the kronecker product of two vectors.

Return the diagonal element of a matrix as a vector

Returns a triangle matrix whose lower diagnonal (including the diagonal) filled with 1. N,M specifies the shape of the return array. K will adjust the sub-diagonal -- positive K fills more 1s.

Returns the copy of matrix with elements above the k-th diagonal zeroed. Positive K fills less 0s.

Returns the copy of matrix with elements below the k-th diagonal zeroed. Positive K fills more 0s.

Returns a matrix whose k-th diagnonal filled with 1. N,M specifies the shape of the return array. K will adjust the sub-diagonal -- positive K moves it upward.

Returns a matrix where M[i,j] == V[i]^(N-j) when increasing is false (default), and M[i,j] == V[i]^j when increasing is true.

Takes an arbitrary number of matrices and performs multiplication, optimizing the order using the shape information and the associativity to minimize the size of the intermediate matrix.

For binary functions

Square function f(x)= x * x.

Logarithmic function. Unlike Common Lisp LOG, it assumes that the input is a positive real, in favor of efficiency. Another reason is that most Common Lisp implementations do not have specialized arrays for positive-only / negative-only floats, thus the runtime type inference cannot deduce whether the resulting array element type have to be a complex or a real.

Logarithmic function. This function behaves exactly like Common Lisp LOG. Unless it can be proven from the input array element type, such as the input being a specialized array of unsigned-int, it takes a conservative approach, i.e., it returns arrays specialized to the complex numbers.

Logarithm of base 2.

This code extends alexandria:shuffle. It additionally accepts arrays and shuffles the elements according to the first axis, viewing the remaining axes as one "element".

Original documentation:

Returns a random permutation of SEQUENCE bounded by START and END. Original sequence may be destructively modified, and (if it contains CONS or lists themselv) share storage with the original one. Signals an error if SEQUENCE is not a proper sequence.

Returns a bit array whose elements are 1 with probability P

pvals is a sequence of probabilities summing up to 1.

Returns a fixnum vector representing a histogram of values. The interval between LOW and HIGH are split by STEP value. All values less than LOW are put in the 0-th bucket; All values greater than equal to HIGH are put in the last bucket.

Takes a multidimentional array and a row-major-index. Returns a list containing the normal index.

Returns a list of the multidimentional indices of the elements which satisfies the predicate FN. Note that the list elements are the multidimentional indices, even for a single-dimensional array.

Returns a list of list of indices of the elements which satisfies the predicate FN. The first list contains the indices for the 1st dimension, the second list contains the indices for the 2nd dimension, and so on.

collect multidimentional indices where the element is nonzero

Collect the elements using a list of multidimentional indices (in a format returned by WHERE).