Generic Math

To remedy the situation of non-extensible math operators in Common Lisp, generic-math provides its own versions of almost all of the standard Lisp math functions/macros. (If the idea seems to be general, that's the rule for inclusion.)

generic-math is used as the framework for tensor, error-propogation, quantity, and histogram as well as being implicitly used in all other parts of cl-ana.

Survey of extensions to Lisp's math

• tensor: Allows sequences/nested sequences to be treated as tensors (multi-indexed arrays). Thus things like (+ (list 2 3 4) 5) evaluate to (7 8 9). Also provides tensor-contract which is a generalization of the dot product for vectors & matrix multiplication.
• error-propogation: Automatically propogates gaussian errors during arithmetic via err-num type. Useful for scientific computing. err-nums are created with the +- function, and are printed back with #~(...) notation so they can be read back via reader macros.
• quantity: Allows mathematical treatment of values with units. Comes with lots of common units and physical constants already defined; also provides a framework for defining new units. Keyword symbols are used to denote units so it is easy to define constants (e.g. speed of light -> (* 299792458 (/ :meter :second))). Metric prefixes are provided as functions, e.g. (centi :meter)
• histogram: As covered in the section on histograms, arithmetic with histograms is provided with the +, -, *, /, and protected-/ functions.

Tensors

Tensors as used in this project are not exactly the same thing as tensors from math/physics. There is no necessity to talk about covariant/contravariant indices, etc. However: These can be enforced by the user, and the tensor library tries to make this easy to understand and use. The main idea behind tensors in this library is to provide a way to use nested sequences as if they are multidimensional arrays. Out of the box, Lisp provides multidimensional arrays, but these are not conveniently used or accessed as if they were simply arrays that contain arrays, and I found in my own experience that this kind of nested array structure would be useful quite often. On top of that, I found it useful to not have to care what type of sequence was being nested to emulate a multidimensional array, so I made the library based on nested sequences rather than a specific underlying type like lists or vectors.

There are two ways to create tensor objects:

1. Use Lisp sequences/nested sequences directly, i.e. '((1 2) (3 4)) is a 4x4 tensor.
2. Use make-tensor to automatically create the sequence/nested sequence.

Rectangular tensors (tensors which have the same size in each dimension no matter where) are safe to use with all tensor functions; non-rectangular ones are only safe to use with tensor-ref and tensor-map as the fancier operations only make sense for rectangular tensors.

Since the first choice is straight forward, we'll show an example of make-tensor:

(defparameter *tensor*
  (make-tensor (list 5 5 5)
               :type 'vector
               :initial-element 0d0))

The keyword arguments :type and :initial-element come by default with the arguments given in the example. make-tensor creates tensors with the same sequence type used at every level; :type selects what this sequence type should be. :initial-element is the default fill value for elements of the tensor.

Tensors have all available generic math functions defined for them automatically as long as the tensor sublibrary is loaded prior to defining the functions with defmath. This also goes for any of your own functions defined with defmath; just make sure you load tensor after defining them and they are added automatically. The default behavior for these methods is to apply them element-wise, and treat non-tensor values as a tensor of the same structure as yours but where each element has the same value.

For example:

(+ (list 1 2 3 4) 1)
==> (2 3 4 5)

This can eliminate the need for looping/mapping over sequences quite often. As a more exotic example:

(+ '((1 2) (3 4))
   '(1 2 3))
==> ((2 3) (5 6))

This happens for because

1. tensor-map treats non-sequence values as tensors of arbitrary depth with the same value everywhere, i.e., referencing a non-sequence always returns the object being referenced no matter how many (or no) indices are applied.
2. tensor-map only maps over indices allowed by the shortest sequence argument. For example: (tensor-map #'+ (list 1 2 3) (list 1 2)) ==> (2 4).

These properties mean that tensor-map can be used on any structure of nested sequences.

Contraction

tensor-contract simultaneously contracts any number of tensors along a single index from each tensor; if you're confused, see http://en.wikipedia.org/wiki/Tensor_contraction (and potentially be more confused). Contraction generalizes both vector dot products and matrix multiplication; as such there is a minimal linear algebra sublibrary in cl-ana.

Example:

(defparameter *mata* '((1 2)
                       (3 4)))
(defparameter *matb* '((5 6)
                       (7 8)))
(tensor-contract (list
                   (cons *mata* 1)
                   (cons *matb* 0)))
==> #(#(19 22) #(43 50))

which is just the nested vector form of the matrix multiple of *mata* and *matb*.

One can implement the physicist's idea of contraction by using a metric tensor along with contraction (it doesn't have to be this ugly, just define your own specific functions to make this nicer):

(defparameter *metric-tensor*
  (function->tensor (list 4 4)
                    (lambda (mu nu)
                      (if (= mu nu)
                          (if (zerop mu)
                              1
                              -1)
                          0))))
(defparameter *four-vector* '(4 1 2 3))
(tensor-contract (list
                   (cons *four-vector* 0)
                   (cons (tensor-contract (list
                                            (cons *four-vector* 0)
                                            (cons *metric-tensor* 0)))
                         0)))
==> 2

This example is made trivial by the lorentz sublibrary, but one special note is the function->tensor function. It takes a function and a list of sizes (along with keyword optional arguments) and generates a nested sequence having the value of the function at each of the possible index-list values.

TODO

How awesome would it be to have algebraic manipulation? And once you define algebraic types/operations you automatically get to use the tensor library with them. I don't have much experience with writing computer algebra systems, but the guys over at maxima do, so I may try to borrow some of their code (it's GPL too).

Next: Dependency Oriented Programming