File	Date	Author	Commit
doc	2015-05-02	Akshay Srinivasan	[11de0a] saving state
examples	2015-12-13	Akshay Srinivasan	[08ec60] Moved field type into MOP. Reorganized class st...
src	2016-08-22	Akshay Srinivasan	[1a015d] update pair function
tests	2016-05-21	Akshay Srinivasan	[6f6003] more bug fixes\!
.gitignore	2011-03-04	Raymond Toy	[ac01ce] Ignore *~.
AUTHORS	2016-03-08	Akshay Srinivasan	[ba99e6] Changed license to Lisp-LGPL
LICENSE	2016-03-08	Akshay Srinivasan	[ba99e6] Changed license to Lisp-LGPL
PREAMBLE	2016-03-08	Akshay Srinivasan	[ba99e6] Changed license to Lisp-LGPL
README-kn.org	2015-12-13	Akshay Srinivasan	[0916da] update kannada doc
README.org	2016-02-14	Akshay Srinivasan	[3688e3] update doc
matlisp.asd	2016-06-19	Akshay Srinivasan	[385d36] added some functions related to polynomials

Read Me

# -*- Mode: org -*-

Common Lisp is a wonderful language for explorative computing, but often lacks tools for numerical computing. This project is an attempt at plugging this perceived shortcoming.

* Overview
  Matlisp provides an infix reader for evaluating mathematical expressions. We use '/' as both a binary
  and unary operator within matlisp's infix reader, this allows us to ensure semantic uniformity even when
  dealing with non-commutativity rings.

  Solving $A x = b$ is consequently as simple as,
  #+BEGIN_SRC lisp
  MATLISP> #i(/A * b)
  #+END_SRC
  , this operation being generic to A being a scalar/matrix/permutation. It's also extremely simple to define
  custom syntax by tweaking "src/reader/infix.lisp".

  For a variety of reasons, Matlisp duplicates general BLAS functionality using Lisp code. These are similar to
  the naive-loop versions found in the BLAS reference implementation, but are generated using macros, and are
  consequently quite terse. These macro expand into lean code blocks which are often very competitive (within
  10-20%) of the corresponding C code. This allows users greater flexibility when dealing with specific
  problems. For instance GEMM is essentially generated by the following code-snippet (exercise: why 'j k i' ?),
  #+BEGIN_SRC lisp
   > (defun mm (A B C)
       (einstein-sum #.(tensor 'double-float) (j k i) (ref C i j) (* (ref A i j) (ref B j k))))
  #+END_SRC
  Matlisp also switches dynamically between BLAS routines (when available) and native lisp
  code to avoid FFI overheads. This functionality can be used (with minimal tweaking) without the need for
  a a BLAS library, or for commutative-ring types which don't have BLAS routines available.

** CLOS, Metaobjects
   Most classes/methods in Matlisp are dynamically generated for specialized types at run-time. Generating class specialized for certain field-types is not particularly interesting, and takes place along usual lines.
   #+BEGIN_SRC lisp
  M> (tensor 'double-float 'dense-tensor)
  |<BLAS-MIXIN DENSE-TENSOR: DOUBLE-FLOAT>|
  M> (tensor 'double-float 'graph-tensor)
  |<GRAPH-TENSOR: DOUBLE-FLOAT>|
   #+END_SRC

   Tying method generation with CLOS, however, comes with its own set of complications. Consider this not-so-hypothetical situation: we're given an abstract class /dense-tensor/, and we wish to write a method for computing the square of a matrix. This sounds simple enough, and as a first pass, one ends up with something like,
  #+BEGIN_SRC lisp
  (defgeneric sqmm (x)
    (:method :before ((x tensor))
       (assert (typep x 'tensor-square-matrix) nil "X must be square")))
  #+END_SRC
  A first solution can be as simple as,
  #+BEGIN_SRC lisp
  (defmethod sqmm ((x dense-tensor) &aux (cl-x (class-name (class-of x))))
    (compile-and-eval
      `(defmethod sqmm ((x ,cl-x) &aux (ret (zeros (dimensions x) (class-of x))))
	 (einstein-sum ,cl-x (j k i) (ref ret i j) (* (ref x i j) (ref x j k)))))
    (sqmm x))
  #+END_SRC
  This works because when a subclass of /dense-tensor/ is encountered which isn't specialized, a method is generated and this method is subsequently called. The next time /sqmm/ is called with an instance of this class, /closer-mop:compute-applicable-methods-using-classes/ places this method at the top of the list. This scheme works perfectly if /dense-tensor/ only has descendents at distance 1 and no more (why ?). Fixing this within the confines of CLOS is troublesome, because every method then needs to check if it is 'special enough'. Cue: Metaobjects to the rescue.

  This is solved in Matlisp in three steps,
  - Define a subclass /tensor-method-generator/ of /closer-mop:funcallable-standard-class/ (a subclass of a function).
  - Define subclasses of /closer-mop:specializer/, one each for marking generated function arguments, and generator arguments.
  - Define custom method dispatch /closer-mop:compute-applicable-methods-using-classes/ for this subclass of generic-functions.
  These definitions can be found in 'src;base;generator.lisp'.

  The solution to the question posed is then simply,
  #+BEGIN_SRC lisp
  (defgeneric sqmm (x)
    (:method :before ((x tensor))
       (assert (typep x 'tensor-square-matrix) nil "X must be square"))
    (:generic-function-class tensor-method-generator))
  (defmethod sqmm ((x #.(make-instance 'group-specializer :object-class 'dense-tensor :group-name :x)) &aux (cl-x (class-name (class-of x))))
    (compile-and-eval
      `(defmethod sqmm ((x ,(make-instance 'classp-specializer :object-class (find-class cl-x))) &aux (ret (zeros (dimensions x) ',cl-x)))
	 (einstein-sum ,cl-x (j k i) (ref ret i j) (* (ref x i j) (ref x j k)))))
    (sqmm x))
  #+END_SRC
  Methods with classp-specializer are only selected by /closer-mop:compute-applicable-methods-using-classes/ if the argument class match that of the classp-specializer instances which are passed as specializers of the method. Similarly, the group-specializer is only selected if all arguments are subclasses of /object-class/, and all argument classes for a given /group-name/ are the same. The CLOS discriminator also ensures that the latter method has very low precedence. This solves the problem posed earlier and code-generating 'Templates' essentially end up becoming first-class functions.

For debuggability and other reasons, this mechanism is abstracted away inside the macro /matlisp:define-tensor-method/.

* Dependencies
- CFFI

  Matlisp uses CFFI for interfacing with foreign libraries. CFFI support varies between lisp implementations;
  Matlisp assumes that certain simple-vectors in the lisp implementation are in unboxed form, and have a SAP
  that can be passed to Fortran functions. It might be possible to replace this with static/pinned arrays, as
  is kind-of done a the moment for foreign dense-tensor.

  Most of the BLAS functionality is duplicated in native lisp code, and so one can very possibly run Matlisp
  without any foreign libraries (but at the cost of losing all LAPACK functionality).

  CFFI is provided by quicklisp.

- closer-MOP

  Matlisp now uses the Meta-Object protocol to generate and maintain classes & methods. This is *not* part of the
  ANSI standard, and consequently only seems to work on SBCL at the moment. Your luck with other implementations
  may wary.

- Trivia

  You may also want to use the latest version of [[https://github.com/guicho271828/trivia/][Trivia]] pattern matching library.
  Matlisp uses the λlist matcher which is a recent addition to the project.

* Install
  Matlisp has been successfully built on SBCL. Getting it running is thematic of quicklisp,
- Install quicklisp http://www.quicklisp.org/beta/.
- Download matlisp source:
#+BEGIN_SRC shell
   > git clone git@github.com:matlisp/matlisp.git
   > ln -s $PWD/matlisp <quicklisp-directory>/local-projects
#+END_SRC
Fire up your lisp implementation and load as usual with quicklisp:
#+BEGIN_SRC lisp
  CL-USER> (ql:quickload :matlisp)
  CL-USER> (in-package :matlisp)
  M> (named-readtables:in-readtable :infix-dispatch-table)
#+END_SRC

* Examples
  #+BEGIN_SRC lisp
  ;;Creation
  M> (copy! (randn '(2 2)) (zeros '(2 2) '((complex double-float))))
  #<|(COMPLEX DOUBLE-FLOAT) STRIDE-ACCESSOR SIMPLE-ARRAY| #(2 2)
    0.8492   -1.976
    2.207    -1.251
  >
  ;;gemv
  M> (let ((a (randn '(2 2)))
	   (b (randn 2)))
       #i(a * b))
  #<|DOUBLE-FLOAT STRIDE-ACCESSOR SIMPLE-ARRAY| #(2)
  1.1885     0.95746
  >

  ;;Tensor contraction
  M> (let ((H (randn '(2 2 2)))
	   (b (randn 2))
	   (c (randn 2))
	   (f (zeros 2)))
	   ;;#i(H @ b @ c)
       (einstein-sum #.(tensor 'double-float) (i j k) (ref f i) (* (ref H i j k) (ref b j) (ref c k))))
  #<|DOUBLE-FLOAT STRIDE-ACCESSOR SIMPLE-ARRAY| #(2)
  0.62586     -1.1128
  >
  ;;Similarly
  M> (let ((H (randn '(2 2 2))))
       #i(H @ randn(2) @ randn(2)))
  #<|DOUBLE-FLOAT STRIDE-ACCESSOR SIMPLE-ARRAY| #(2)
   0.3234  -0.6201
  >
  #+END_SRC

* Documentation
  More documentation will be added as the project reaches a stable state.

* Enchancements
- [[https://github.com/matlisp/matlisp-forbi][matlisp-forbi]]

  The API for BLAS functions dot ensures inconsistent ABIs between compilers. This package provides a Fortran wrapper (and Lisp methods for `dot`) that fixes these issues. It also provides F77 methods for elementwise division, which follow the `scal` API.e

- Weyl

  Weyl is a CAS written in Lisp (and for Lisp!) at Cornell by Richard Zippel's group. Currently, this used only within
  'src;base;symbolic.lisp' (and assoc. infix readers), for working with symbolic expressions. In order to use this functionality,
  Weyl must be loaded before Matlisp.

  [[https://github.com/matlisp/weyl][Weyl]] can installed from 'git@github.com:matlisp/weyl.git'.

* Tracker
** Completed
   * Dynamic class/method generation, using MOP
   * Complete BLAS/LAPACK functionality for types double-float, single-float, (complex single-float), (complex double-float).
   * Partial support for dense-tensor with a foreign-pointer store.
   * Inplace slicing, real-imaginary part views, negative strides for dense-tensors.
   * permutations, sorting, conversion between action/cycle/flip representations.
   * Optimized loop generators (einstein/iterate for-mod) in Lisp; BLAS functionality duplicated, and switches automatically b/w Lisp and Fortran.
   * Arbitrary tensor contraction.
   * Graphs: a general CCS/CCR matrix implementation, lisp adjacency list support, iterate macros for DFS/BFS/SFD graph-traversal, tree-decomposition,
     cholesky-covers, maximum acyclic subgraph, Djikstra's algorithms.
   * Data structures: Fibonacci heap, Union-Find structure, minimal Doubly linked lists.
   * Hash-table sparse tensor: O(1) read/write.
   * Co-ordinate sparse tensor

** TODO Incomplete/Planned
- Random distributions
  Implementing a fast, automated version of [[http://www.jstatsoft.org/article/view/v005i08/ziggurat.pdf][Ziggurat]] algorithm. Given a sampler for the tail of the distribution and the
  form of the density function, it should be theoretically possible to generate a Ziggurat sampler.
- Unify slicing syntax
  Unify the slicing syntax used by iterate for-mod/einstein macros. Unify these with a more powerful language.
- Automatic Differentiation
- Symbolic Integration
  Needs extensive hacking of Weyl.
- Gnuplot interface
- (C)Python-bridge
- FiveAM tests

* Emacs
Matlisp uses a variety of Unicode symbols for some function names and certain operators in the infix reader.
The user can readily change these to his suiting, or instead use the following Emacs shortcuts to enter these
characters.

#+BEGIN_SRC lisp
;; Lisp
(defun add-lisp-slime-hook (func)
  (add-hook 'lisp-mode-hook func)
  (add-hook 'slime-repl-mode-hook func))
;;#\GREEK_SMALL_LETTER_LAMDA is bound to lambda in :infix-dispatch-table; inherited from λ-reader
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c \\") (lambda () (interactive (insert "λ"))))))
;;#\LONG_RIGHTWARDS_ARROW_FROM_BAR used for anonymous function definition in Infix
;;#i([x] ⟼ x + 1)
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c /") (lambda () (interactive (insert "⟼"))))))
;;#\CIRCLE_TIMES used for tensor-product in Infix
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c *") (lambda () (interactive (insert "⊗"))))))
;;#\MIDDLE_DOT used for tensor-contraction (also bound is @) in Infix
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c .") (lambda () (interactive (insert "·"))))))
;;Used in the function `(δ-i g &optional i j)` in graph-accessor.lisp
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c a") (lambda () (interactive (insert "δ"))))))
;;#\DEVANAGARI_LETTER_SA used in infix for tensors involving symbolic expressions.
(add-lisp-slime-hook #'(lambda () (local-set-key (kbd "C-c s") (lambda () (interactive (insert "स"))))))
#+END_SRC