1099 lines
38 KiB
HTML
1099 lines
38 KiB
HTML
|
<!DOCTYPE html>
|
|||
|
|
<html lang="en">
|
|||
|
|
<head>
|
|||
|
|
<meta name="generator" content=
|
|||
|
|
"HTML Tidy for HTML5 for Linux version 5.2.0">
|
|||
|
|
<title>Multidimensional arrays</title>
|
|||
|
|
<meta charset="utf-8">
|
|||
|
|
<meta name="description" content="A collection of examples of using Common Lisp">
|
|||
|
|
<meta name="viewport" content=
|
|||
|
|
"width=device-width, initial-scale=1">
|
|||
|
|
<link rel="stylesheet" href=
|
|||
|
|
"assets/style.css">
|
|||
|
|
<script type="text/javascript" src=
|
|||
|
|
"assets/highlight-lisp.js">
|
|||
|
|
</script>
|
|||
|
|
<script type="text/javascript" src=
|
|||
|
|
"assets/jquery-3.2.1.min.js">
|
|||
|
|
</script>
|
|||
|
|
<script type="text/javascript" src=
|
|||
|
|
"assets/jquery.toc/jquery.toc.min.js">
|
|||
|
|
</script>
|
|||
|
|
<script type="text/javascript" src=
|
|||
|
|
"assets/toggle-toc.js">
|
|||
|
|
</script>
|
|||
|
|
|
|||
|
|
<link rel="stylesheet" href=
|
|||
|
|
"assets/github.css">
|
|||
|
|
|
|||
|
|
<link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css" integrity="sha384-BVYiiSIFeK1dGmJRAkycuHAHRg32OmUcww7on3RYdg4Va+PmSTsz/K68vbdEjh4u" crossorigin="anonymous">
|
|||
|
|
</head>
|
|||
|
|
<body>
|
|||
|
|
<h1 id="title-xs"><a href="index.html">The Common Lisp Cookbook</a> – Multidimensional arrays</h1>
|
|||
|
|
<div id="logo-container">
|
|||
|
|
<a href="index.html">
|
|||
|
|
<img id="logo" src="assets/cl-logo-blue.png"/>
|
|||
|
|
</a>
|
|||
|
|
|
|||
|
|
<div id="searchform-container">
|
|||
|
|
<form onsubmit="duckSearch()" action="javascript:void(0)">
|
|||
|
|
<input id="searchField" type="text" value="" placeholder="Search...">
|
|||
|
|
</form>
|
|||
|
|
</div>
|
|||
|
|
|
|||
|
|
<div id="toc-container" class="toc-close">
|
|||
|
|
<div id="toc-title">Table of Contents</div>
|
|||
|
|
<ul id="toc" class="list-unstyled"></ul>
|
|||
|
|
</div>
|
|||
|
|
</div>
|
|||
|
|
|
|||
|
|
<div id="content-container">
|
|||
|
|
<h1 id="title-non-xs"><a href="index.html">The Common Lisp Cookbook</a> – Multidimensional arrays</h1>
|
|||
|
|
|
|||
|
|
<!-- Announcement we can keep for 1 month or more. I remove it and re-add it from time to time. -->
|
|||
|
|
<p class="announce">
|
|||
|
|
📹 <a href="https://www.udemy.com/course/common-lisp-programming/?couponCode=6926D599AA-LISP4ALL">NEW! Learn Lisp in videos and support our contributors with this 40% discount.</a>
|
|||
|
|
</p>
|
|||
|
|
<p class="announce-neutral">
|
|||
|
|
📕 <a href="index.html#download-in-epub">Get the EPUB and PDF</a>
|
|||
|
|
</p>
|
|||
|
|
|
|||
|
|
|
|||
|
|
<div id="content"
|
|||
|
|
<p>Common Lisp has native support for multidimensional arrays, with some
|
|||
|
|
special treatment for 1-D arrays, called <code>vectors</code>. Arrays can be
|
|||
|
|
<em>generalised</em> and contain any type (<code>element-type t</code>), or they
|
|||
|
|
can be <em>specialised</em> to contain specific types such as <code>single-float</code>
|
|||
|
|
or <code>integer</code>. A good place to start is
|
|||
|
|
<a href="http://www.gigamonkeys.com/book/collections.html">Practical Common Lisp Chapter 11, Collections</a> by
|
|||
|
|
Peter Seibel.</p>
|
|||
|
|
|
|||
|
|
<p>A quick reference to some common operations on arrays is given in the
|
|||
|
|
section on <a href="data-structures.html">Arrays and vectors</a>.</p>
|
|||
|
|
|
|||
|
|
<p>Some libraries available on <a href="https://www.quicklisp.org/beta/">Quicklisp</a> for manipulating arrays:</p>
|
|||
|
|
|
|||
|
|
<ul>
|
|||
|
|
<li><a href="https://github.com/tpapp/array-operations">array-operations</a> defines
|
|||
|
|
functions <code>generate</code>, <code>permute</code>, <code>displace</code>, <code>flatten</code>, <code>split</code>,
|
|||
|
|
<code>combine</code>, <code>reshape</code>. It also defines <code>each</code>, for element-wise
|
|||
|
|
operations. This library is not maintained by the original author,
|
|||
|
|
but there is an <a href="https://github.com/bendudson/array-operations">actively maintained fork</a>.</li>
|
|||
|
|
<li><a href="https://github.com/rigetticomputing/cmu-infix">cmu-infix</a> includes
|
|||
|
|
array indexing syntax for multidimensional arrays.</li>
|
|||
|
|
<li><a href="https://github.com/tpapp/lla">lla</a> is a library for linear algebra, calling BLAS and LAPACK
|
|||
|
|
libraries. It differs from most CL linear algebra packages in using
|
|||
|
|
intuitive function names, and can operate on native arrays as well as
|
|||
|
|
CLOS objects.</li>
|
|||
|
|
</ul>
|
|||
|
|
|
|||
|
|
<p>This page covers what can be done with the built-in multidimensional
|
|||
|
|
arrays, but there are limitations. In particular:</p>
|
|||
|
|
|
|||
|
|
<ul>
|
|||
|
|
<li>Interoperability with foreign language arrays, for example when
|
|||
|
|
calling libraries such as BLAS, LAPACK or GSL.</li>
|
|||
|
|
<li>Extending arithmetic and other mathematical operators to handle
|
|||
|
|
arrays, for example so that <code>(+ a b)</code> works
|
|||
|
|
when <code>a</code> and/or <code>b</code> are arrays.</li>
|
|||
|
|
</ul>
|
|||
|
|
|
|||
|
|
<p>Both of these problems can be solved by using CLOS to define an
|
|||
|
|
extended array class, with native arrays as a special case.
|
|||
|
|
Some libraries available through
|
|||
|
|
<a href="https://www.quicklisp.org/beta/">quicklisp</a> which take this approach
|
|||
|
|
are:</p>
|
|||
|
|
|
|||
|
|
<ul>
|
|||
|
|
<li><a href="https://github.com/bharath1097/matlisp/">matlisp</a>, some of which is
|
|||
|
|
described in sections below.</li>
|
|||
|
|
<li><a href="https://github.com/melisgl/mgl-mat">MGL-MAT</a>, which has a manual
|
|||
|
|
and provides bindings to BLAS and CUDA. This is used in a machine
|
|||
|
|
learning library <a href="https://github.com/melisgl/mgl">MGL</a>.</li>
|
|||
|
|
<li><a href="https://github.com/ghollisjr/cl-ana/wiki">cl-ana</a>, a data analysis
|
|||
|
|
package with a manual, which includes operations on arrays.</li>
|
|||
|
|
<li><a href="https://www.common-lisp.net/project/antik/">Antik</a>, used in
|
|||
|
|
<a href="https://common-lisp.net/project/gsll/">GSLL</a>, a binding to the GNU
|
|||
|
|
Scientific Library.</li>
|
|||
|
|
</ul>
|
|||
|
|
|
|||
|
|
<p>A relatively new but actively developed package is
|
|||
|
|
<a href="https://github.com/rigetticomputing/magicl">MAGICL</a>, which provides
|
|||
|
|
wrappers around BLAS and LAPACK libraries. At the time of writing this
|
|||
|
|
package is not on Quicklisp, and only works under SBCL and CCL. It
|
|||
|
|
seems to be particularly focused on complex arrays, but not
|
|||
|
|
exclusively.
|
|||
|
|
To install, clone the repository in your quicklisp <code>local-projects</code>
|
|||
|
|
directory e.g. under Linux/Unix:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-bash">$ cd ~/quicklisp/local-projects
|
|||
|
|
$ git clone https://github.com/rigetticomputing/magicl.git
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Instructions for installing dependencies (BLAS, LAPACK and Expokit)
|
|||
|
|
are given on the <a href="https://github.com/rigetticomputing/magicl">github web pages</a>.
|
|||
|
|
Low-level routines wrap foreign functions, so have the Fortran names
|
|||
|
|
e.g <code>magicl.lapack-cffi::%zgetrf</code>. Higher-level interfaces to some of
|
|||
|
|
these functions also exist, see the
|
|||
|
|
<a href="https://github.com/rigetti/magicl/blob/master/src/high-level/">source directory</a> and <a href="https://github.com/quil-lang/magicl/blob/master/doc/high-level.md">documentation</a>.</p>
|
|||
|
|
|
|||
|
|
<p>Taking this further, domain specific languages have been built on Common
|
|||
|
|
Lisp, which can be used for numerical calculations with arrays.
|
|||
|
|
At the time of writing the most widely used and supported of these are:</p>
|
|||
|
|
|
|||
|
|
<ul>
|
|||
|
|
<li><a href="http://maxima.sourceforge.net/documentation.html">Maxima</a></li>
|
|||
|
|
<li><a href="https://github.com/daly/axiom">Axiom</a></li>
|
|||
|
|
</ul>
|
|||
|
|
|
|||
|
|
<p><a href="https://github.com/drmeister/clasp">CLASP</a> is a project
|
|||
|
|
which aims to ease interoperability of Common Lisp with other
|
|||
|
|
languages (particularly C++), by using <a href="http://llvm.org/">LLVM</a>.
|
|||
|
|
One of the main applications of this project is to numerical/scientific
|
|||
|
|
computing.</p>
|
|||
|
|
|
|||
|
|
<h2 id="creating">Creating</h2>
|
|||
|
|
|
|||
|
|
<p>The function <a href="http://clhs.lisp.se/Body/f_mk_ar.htm">CLHS: make-array</a>
|
|||
|
|
can create arrays filled with a single value</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter *my-array* (make-array '(3 2) :initial-element 1.0))
|
|||
|
|
*MY-ARRAY*
|
|||
|
|
* *my-array*
|
|||
|
|
#2A((1.0 1.0) (1.0 1.0) (1.0 1.0))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>More complicated array values can be generated by first making an
|
|||
|
|
array, and then iterating over the elements to fill in the values (see
|
|||
|
|
section below on element access).</p>
|
|||
|
|
|
|||
|
|
<p>The <a href="https://github.com/tpapp/array-operations">array-operations</a>
|
|||
|
|
library provides <code>generate</code>, a convenient
|
|||
|
|
function for creating arrays which wraps this iteration.</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (ql:quickload :array-operations)
|
|||
|
|
To load "array-operations":
|
|||
|
|
Load 1 ASDF system:
|
|||
|
|
array-operations
|
|||
|
|
; Loading "array-operations"
|
|||
|
|
|
|||
|
|
(:ARRAY-OPERATIONS)
|
|||
|
|
|
|||
|
|
* (aops:generate #'identity 7 :position)
|
|||
|
|
#(0 1 2 3 4 5 6)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Note that the nickname for <code>array-operations</code> is <code>aops</code>. The
|
|||
|
|
<code>generate</code> function can also iterate over the array subscripts by
|
|||
|
|
passing the key <code>:subscripts</code>. See the
|
|||
|
|
<a href="https://github.com/tpapp/array-operations">github repository</a> for
|
|||
|
|
more examples.</p>
|
|||
|
|
|
|||
|
|
<h3 id="random-numbers">Random numbers</h3>
|
|||
|
|
|
|||
|
|
<p>To create an 3x3 array containing random numbers drawn from a uniform
|
|||
|
|
distribution, <code>generate</code> can be used to call the CL
|
|||
|
|
<a href="http://clhs.lisp.se/Body/f_random.htm">random</a> function:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (aops:generate (lambda () (random 1.0)) '(3 3))
|
|||
|
|
#2A((0.99292254 0.929777 0.93538976)
|
|||
|
|
(0.31522608 0.45167792 0.9411855)
|
|||
|
|
(0.96221936 0.9143338 0.21972346))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>An array of Gaussian (normal) random numbers with mean of zero and
|
|||
|
|
standard deviation of one, using the
|
|||
|
|
<a href="https://common-lisp.net/project/alexandria/">alexandria</a> package:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (ql:quickload :alexandria)
|
|||
|
|
To load "alexandria":
|
|||
|
|
Load 1 ASDF system:
|
|||
|
|
alexandria
|
|||
|
|
; Loading "alexandria"
|
|||
|
|
|
|||
|
|
(:ALEXANDRIA)
|
|||
|
|
|
|||
|
|
* (aops:generate #'alexandria:gaussian-random 4)
|
|||
|
|
#(0.5522547885338768d0 -1.2564808468164517d0 0.9488161476129733d0
|
|||
|
|
-0.10372852118266523d0)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Note that this is not particularly efficient: It requires a function
|
|||
|
|
call for each element, and although <code>gaussian-random</code> returns
|
|||
|
|
two random numbers, only one of them is used.</p>
|
|||
|
|
|
|||
|
|
<p>For more efficient implementations, and a wider range of probability
|
|||
|
|
distributions, there are packages available on Quicklisp. See
|
|||
|
|
<a href="https://www.cliki.net/statistics">CLiki</a> for a list.</p>
|
|||
|
|
|
|||
|
|
<h2 id="accessing-elements">Accessing elements</h2>
|
|||
|
|
|
|||
|
|
<p>To access the individual elements of an array there are the <a href="http://clhs.lisp.se/Body/f_aref.htm">aref</a>
|
|||
|
|
and <a href="http://clhs.lisp.se/Body/f_row_ma.htm#row-major-aref">row-major-aref</a> functions.</p>
|
|||
|
|
|
|||
|
|
<p>The <a href="http://clhs.lisp.se/Body/f_aref.htm">aref</a> function takes the same number of index arguments as the
|
|||
|
|
array has dimensions. Indexing is from 0 and row-major as in C, but
|
|||
|
|
not Fortran.</p>
|
|||
|
|
<pre><code class="language-lisp">* (defparameter *a* #(1 2 3 4))
|
|||
|
|
*A*
|
|||
|
|
* (aref *a* 0)
|
|||
|
|
1
|
|||
|
|
* (aref *a* 3)
|
|||
|
|
4
|
|||
|
|
* (defparameter *b* #2A((1 2 3) (4 5 6)))
|
|||
|
|
*B*
|
|||
|
|
* (aref *b* 1 0)
|
|||
|
|
4
|
|||
|
|
* (aref *b* 0 2)
|
|||
|
|
3
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>The range of these indices can be found using
|
|||
|
|
<a href="http://clhs.lisp.se/Body/f_ar_d_1.htm">array-dimensions</a>:</p>
|
|||
|
|
|
|||
|
|
<pre><code>* (array-dimensions *a*)
|
|||
|
|
(4)
|
|||
|
|
* (array-dimensions *b*)
|
|||
|
|
(2 3)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>or the rank of the array can be found, and then the size of each
|
|||
|
|
dimension queried:</p>
|
|||
|
|
<pre><code class="language-lisp">* (array-rank *a*)
|
|||
|
|
1
|
|||
|
|
* (array-dimension *a* 0)
|
|||
|
|
4
|
|||
|
|
* (array-rank *b*)
|
|||
|
|
2
|
|||
|
|
* (array-dimension *b* 0)
|
|||
|
|
2
|
|||
|
|
* (array-dimension *b* 1)
|
|||
|
|
3
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>To loop over an array nested loops can be used, such as</p>
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a #2A((1 2 3) (4 5 6)))
|
|||
|
|
A
|
|||
|
|
* (destructuring-bind (n m) (array-dimensions a)
|
|||
|
|
(loop for i from 0 below n do
|
|||
|
|
(loop for j from 0 below m do
|
|||
|
|
(format t "a[~a ~a] = ~a~%" i j (aref a i j)))))
|
|||
|
|
|
|||
|
|
a[0 0] = 1
|
|||
|
|
a[0 1] = 2
|
|||
|
|
a[0 2] = 3
|
|||
|
|
a[1 0] = 4
|
|||
|
|
a[1 1] = 5
|
|||
|
|
a[1 2] = 6
|
|||
|
|
NIL
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>A utility macro which does this for multiple dimensions is <code>nested-loop</code>:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">(defmacro nested-loop (syms dimensions &body body)
|
|||
|
|
"Iterates over a multidimensional range of indices.
|
|||
|
|
|
|||
|
|
SYMS must be a list of symbols, with the first symbol
|
|||
|
|
corresponding to the outermost loop.
|
|||
|
|
|
|||
|
|
DIMENSIONS will be evaluated, and must be a list of
|
|||
|
|
dimension sizes, of the same length as SYMS.
|
|||
|
|
|
|||
|
|
Example:
|
|||
|
|
(nested-loop (i j) '(10 20) (format t '~a ~a~%' i j))
|
|||
|
|
|
|||
|
|
"
|
|||
|
|
(unless syms (return-from nested-loop `(progn ,@body))) ; No symbols
|
|||
|
|
|
|||
|
|
;; Generate gensyms for dimension sizes
|
|||
|
|
(let* ((rank (length syms))
|
|||
|
|
(syms-rev (reverse syms)) ; Reverse, since starting with innermost
|
|||
|
|
(dims-rev (loop for i from 0 below rank collecting (gensym))) ; innermost dimension first
|
|||
|
|
(result `(progn ,@body))) ; Start with innermost expression
|
|||
|
|
;; Wrap previous result inside a loop for each dimension
|
|||
|
|
(loop for sym in syms-rev for dim in dims-rev do
|
|||
|
|
(unless (symbolp sym) (error "~S is not a symbol. First argument to nested-loop must be a list of symbols" sym))
|
|||
|
|
(setf result
|
|||
|
|
`(loop for ,sym from 0 below ,dim do
|
|||
|
|
,result)))
|
|||
|
|
;; Add checking of rank and dimension types, and get dimensions into gensym list
|
|||
|
|
(let ((dims (gensym)))
|
|||
|
|
`(let ((,dims ,dimensions))
|
|||
|
|
(unless (= (length ,dims) ,rank) (error "Incorrect number of dimensions: Expected ~a but got ~a" ,rank (length ,dims)))
|
|||
|
|
(dolist (dim ,dims)
|
|||
|
|
(unless (integerp dim) (error "Dimensions must be integers: ~S" dim)))
|
|||
|
|
(destructuring-bind ,(reverse dims-rev) ,dims ; Dimensions reversed so that innermost is last
|
|||
|
|
,result)))))
|
|||
|
|
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>so that the contents of a 2D array can be printed using</p>
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a #2A((1 2 3) (4 5 6)))
|
|||
|
|
A
|
|||
|
|
* (nested-loop (i j) (array-dimensions a)
|
|||
|
|
(format t "a[~a ~a] = ~a~%" i j (aref a i j)))
|
|||
|
|
|
|||
|
|
a[0 0] = 1
|
|||
|
|
a[0 1] = 2
|
|||
|
|
a[0 2] = 3
|
|||
|
|
a[1 0] = 4
|
|||
|
|
a[1 1] = 5
|
|||
|
|
a[1 2] = 6
|
|||
|
|
NIL
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>[Note: This macro is available in <a href="https://github.com/bendudson/array-operations">this fork</a> of array-operations, but
|
|||
|
|
not Quicklisp]</p>
|
|||
|
|
|
|||
|
|
<h3 id="row-major-indexing">Row major indexing</h3>
|
|||
|
|
|
|||
|
|
<p>In some cases, particularly element-wise operations, the number of
|
|||
|
|
dimensions does not matter. To write code which is independent of the
|
|||
|
|
number of dimensions, array element access can be done using a
|
|||
|
|
single flattened index via
|
|||
|
|
<a href="http://clhs.lisp.se/Body/f_row_ma.htm#row-major-aref">row-major-aref</a>.
|
|||
|
|
The array size is given by <a href="http://clhs.lisp.se/Body/f_ar_tot.htm">array-total-size</a>, with the flattened
|
|||
|
|
index starting at 0.</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a #2A((1 2 3) (4 5 6)))
|
|||
|
|
A
|
|||
|
|
* (array-total-size a)
|
|||
|
|
6
|
|||
|
|
* (loop for i from 0 below (array-total-size a) do
|
|||
|
|
(setf (row-major-aref a i) (+ 2.0 (row-major-aref a i))))
|
|||
|
|
NIL
|
|||
|
|
* a
|
|||
|
|
#2A((3.0 4.0 5.0) (6.0 7.0 8.0))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h3 id="infix-syntax">Infix syntax</h3>
|
|||
|
|
|
|||
|
|
<p>The <a href="https://github.com/rigetticomputing/cmu-infix">cmu-infix</a> library
|
|||
|
|
provides some different syntax which can make mathematical expressions
|
|||
|
|
easier to read:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (ql:quickload :cmu-infix)
|
|||
|
|
To load "cmu-infix":
|
|||
|
|
Load 1 ASDF system:
|
|||
|
|
cmu-infix
|
|||
|
|
; Loading "cmu-infix"
|
|||
|
|
|
|||
|
|
(:CMU-INFIX)
|
|||
|
|
|
|||
|
|
* (named-readtables:in-readtable cmu-infix:syntax)
|
|||
|
|
(("COMMON-LISP-USER" . #<NAMED-READTABLE CMU-INFIX:SYNTAX {10030158B3}>)
|
|||
|
|
...)
|
|||
|
|
|
|||
|
|
* (defparameter arr (make-array '(3 2) :initial-element 1.0))
|
|||
|
|
ARR
|
|||
|
|
|
|||
|
|
* #i(arr[0 1] = 2.0)
|
|||
|
|
2.0
|
|||
|
|
|
|||
|
|
* arr
|
|||
|
|
#2A((1.0 2.0) (1.0 1.0) (1.0 1.0))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>A matrix-matrix multiply operation can be implemented as:</p>
|
|||
|
|
<pre><code class="language-lisp">(let ((A #2A((1 2) (3 4)))
|
|||
|
|
(B #2A((5 6) (7 8)))
|
|||
|
|
(result (make-array '(2 2) :initial-element 0.0)))
|
|||
|
|
|
|||
|
|
(loop for i from 0 to 1 do
|
|||
|
|
(loop for j from 0 to 1 do
|
|||
|
|
(loop for k from 0 to 1 do
|
|||
|
|
#i(result[i j] += A[i k] * B[k j]))))
|
|||
|
|
result)
|
|||
|
|
</code></pre>
|
|||
|
|
<p>See the section below on linear algebra, for alternative
|
|||
|
|
matrix-multiply implementations.</p>
|
|||
|
|
|
|||
|
|
<h2 id="element-wise-operations">Element-wise operations</h2>
|
|||
|
|
|
|||
|
|
<p>To multiply two arrays of numbers of the same size, pass a function
|
|||
|
|
to <code>each</code> in the <a href="https://github.com/tpapp/array-operations">array-operations</a> library:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (aops:each #'* #(1 2 3) #(2 3 4))
|
|||
|
|
#(2 6 12)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>For improved efficiency there is the <code>aops:each*</code> function, which
|
|||
|
|
takes a type as first argument to specialise the result array.</p>
|
|||
|
|
|
|||
|
|
<p>To add a constant to all elements of an array:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter *a* #(1 2 3 4))
|
|||
|
|
*A*
|
|||
|
|
* (aops:each (lambda (it) (+ 42 it)) *a*)
|
|||
|
|
#(43 44 45 46)
|
|||
|
|
* *a*
|
|||
|
|
#(1 2 3 4)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Note that <code>each</code> is not destructive, but makes a new array.
|
|||
|
|
All arguments to <code>each</code> must be arrays of the same size,
|
|||
|
|
so <code>(aops:each #'+ 42 *a*)</code> is not valid.</p>
|
|||
|
|
|
|||
|
|
<h3 id="vectorising-expressions">Vectorising expressions</h3>
|
|||
|
|
|
|||
|
|
<p>An alternative approach to the <code>each</code> function above, is to use a
|
|||
|
|
macro to iterate over all elements of an array:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">(defmacro vectorize (variables &body body)
|
|||
|
|
;; Check that variables is a list of only symbols
|
|||
|
|
(dolist (var variables)
|
|||
|
|
(if (not (symbolp var))
|
|||
|
|
(error "~S is not a symbol" var)))
|
|||
|
|
|
|||
|
|
;; Get the size of the first variable, and create a new array
|
|||
|
|
;; of the same type for the result
|
|||
|
|
`(let ((size (array-total-size ,(first variables))) ; Total array size (same for all variables)
|
|||
|
|
(result (make-array (array-dimensions ,(first variables)) ; Returned array
|
|||
|
|
:element-type (array-element-type ,(first variables)))))
|
|||
|
|
;; Check that all variables have the same sizeo
|
|||
|
|
,@(mapcar (lambda (var) `(if (not (equal (array-dimensions ,(first variables))
|
|||
|
|
(array-dimensions ,var)))
|
|||
|
|
(error "~S and ~S have different dimensions" ',(first variables) ',var)))
|
|||
|
|
(rest variables))
|
|||
|
|
|
|||
|
|
(dotimes (indx size)
|
|||
|
|
;; Locally redefine variables to be scalars at a given index
|
|||
|
|
(let ,(mapcar (lambda (var) (list var `(row-major-aref ,var indx))) variables)
|
|||
|
|
;; User-supplied function body now evaluated for each index in turn
|
|||
|
|
(setf (row-major-aref result indx) (progn ,@body))))
|
|||
|
|
result))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>[Note: Expanded versions of this macro are available in <a href="https://github.com/bendudson/array-operations">this
|
|||
|
|
fork</a> of array-operations, but
|
|||
|
|
not Quicklisp]</p>
|
|||
|
|
|
|||
|
|
<p>This can be used as:</p>
|
|||
|
|
<pre><code class="language-lisp">* (defparameter *a* #(1 2 3 4))
|
|||
|
|
*A*
|
|||
|
|
* (vectorize (*a*) (* 2 *a*))
|
|||
|
|
#(2 4 6 8)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Inside the body of the expression (second form in <code>vectorize</code>
|
|||
|
|
expression) the symbol <code>*a*</code> is bound to a single element.
|
|||
|
|
This means that the built-in mathematical functions can be used:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a #(1 2 3 4))
|
|||
|
|
A
|
|||
|
|
* (defparameter b #(2 3 4 5))
|
|||
|
|
B
|
|||
|
|
* (vectorize (a b) (* a (sin b)))
|
|||
|
|
#(0.9092974 0.28224 -2.2704074 -3.8356972)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>and combined with <code>cmu-infix</code></p>
|
|||
|
|
<pre><code class="language-lisp">* (vectorize (a b) #i(a * sin(b)) )
|
|||
|
|
#(0.9092974 0.28224 -2.2704074 -3.8356972)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h3 id="calling-blas">Calling BLAS</h3>
|
|||
|
|
|
|||
|
|
<p>Several packages provide wrappers around BLAS, for fast matrix manipulation.</p>
|
|||
|
|
|
|||
|
|
<p>The <a href="https://github.com/tpapp/lla">lla</a> package in quicklisp includes
|
|||
|
|
calls to some functions:</p>
|
|||
|
|
|
|||
|
|
<h4 id="scale-an-array">Scale an array</h4>
|
|||
|
|
|
|||
|
|
<p>scaling by a constant factor:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a #(1 2 3))
|
|||
|
|
* (lla:scal! 2.0 a)
|
|||
|
|
* a
|
|||
|
|
#(2.0d0 4.0d0 6.0d0)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h4 id="axpy">AXPY</h4>
|
|||
|
|
|
|||
|
|
<p>This calculates <code>a * x + y</code> where <code>a</code> is a constant, <code>x</code> and <code>y</code> are arrays.
|
|||
|
|
The <code>lla:axpy!</code> function is destructive, modifying the last argument (<code>y</code>).</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter x #(1 2 3))
|
|||
|
|
A
|
|||
|
|
* (defparameter y #(2 3 4))
|
|||
|
|
B
|
|||
|
|
* (lla:axpy! 0.5 x y)
|
|||
|
|
#(2.5d0 4.0d0 5.5d0)
|
|||
|
|
* x
|
|||
|
|
#(1.0d0 2.0d0 3.0d0)
|
|||
|
|
* y
|
|||
|
|
#(2.5d0 4.0d0 5.5d0)
|
|||
|
|
</code></pre>
|
|||
|
|
<p>If the <code>y</code> array is complex, then this operation calls the complex
|
|||
|
|
number versions of these operators:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter x #(1 2 3))
|
|||
|
|
* (defparameter y (make-array 3 :element-type '(complex double-float)
|
|||
|
|
:initial-element #C(1d0 1d0)))
|
|||
|
|
* y
|
|||
|
|
#(#C(1.0d0 1.0d0) #C(1.0d0 1.0d0) #C(1.0d0 1.0d0))
|
|||
|
|
|
|||
|
|
* (lla:axpy! #C(0.5 0.5) a b)
|
|||
|
|
#(#C(1.5d0 1.5d0) #C(2.0d0 2.0d0) #C(2.5d0 2.5d0))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h4 id="dot-product">Dot product</h4>
|
|||
|
|
|
|||
|
|
<p>The dot product of two vectors:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter x #(1 2 3))
|
|||
|
|
* (defparameter y #(2 3 4))
|
|||
|
|
* (lla:dot x y)
|
|||
|
|
20.0d0
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h3 id="reductions">Reductions</h3>
|
|||
|
|
|
|||
|
|
<p>The <a href="http://clhs.lisp.se/Body/f_reduce.htm">reduce</a> function operates
|
|||
|
|
on sequences, including vectors (1D arrays), but not on
|
|||
|
|
multidimensional arrays.
|
|||
|
|
To get around this, multidimensional arrays can be displaced to create
|
|||
|
|
a 1D vector.
|
|||
|
|
Displaced arrays share storage with the original array, so this is a
|
|||
|
|
fast operation which does not require copying data:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a #2A((1 2) (3 4)))
|
|||
|
|
A
|
|||
|
|
* (reduce #'max (make-array (array-total-size a) :displaced-to a))
|
|||
|
|
4
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>The <code>array-operations</code> package contains <code>flatten</code>, which returns a
|
|||
|
|
displaced array i.e doesn’t copy data:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (reduce #'max (aops:flatten a))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>An SBCL extension,
|
|||
|
|
<a href="http://www.sbcl.org/manual/#index-array_002dstorage_002dvector">array-storage-vector</a>
|
|||
|
|
provides an efficient but not portable way to achieve the same thing:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (reduce #'max (array-storage-vector a))
|
|||
|
|
4
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>More complex reductions are sometimes needed, for example finding the
|
|||
|
|
maximum absolute difference between two arrays. Using the above
|
|||
|
|
methods we could do:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a #2A((1 2) (3 4)))
|
|||
|
|
A
|
|||
|
|
* (defparameter b #2A((1 3) (5 4)))
|
|||
|
|
B
|
|||
|
|
* (reduce #'max (aops:flatten
|
|||
|
|
(aops:each
|
|||
|
|
(lambda (a b) (abs (- a b))) a b)))
|
|||
|
|
2
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>This involves allocating an array to hold the intermediate result,
|
|||
|
|
which for large arrays could be inefficient. Similarly to <code>vectorize</code>
|
|||
|
|
defined above, a macro which does not allocate can be defined as:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">(defmacro vectorize-reduce (fn variables &body body)
|
|||
|
|
"Performs a reduction using FN over all elements in a vectorized expression
|
|||
|
|
on array VARIABLES.
|
|||
|
|
|
|||
|
|
VARIABLES must be a list of symbols bound to arrays.
|
|||
|
|
Each array must have the same dimensions. These are
|
|||
|
|
checked at compile and run-time respectively.
|
|||
|
|
"
|
|||
|
|
;; Check that variables is a list of only symbols
|
|||
|
|
(dolist (var variables)
|
|||
|
|
(if (not (symbolp var))
|
|||
|
|
(error "~S is not a symbol" var)))
|
|||
|
|
|
|||
|
|
(let ((size (gensym)) ; Total array size (same for all variables)
|
|||
|
|
(result (gensym)) ; Returned value
|
|||
|
|
(indx (gensym))) ; Index inside loop from 0 to size
|
|||
|
|
|
|||
|
|
;; Get the size of the first variable
|
|||
|
|
`(let ((,size (array-total-size ,(first variables))))
|
|||
|
|
;; Check that all variables have the same size
|
|||
|
|
,@(mapcar (lambda (var) `(if (not (equal (array-dimensions ,(first variables))
|
|||
|
|
(array-dimensions ,var)))
|
|||
|
|
(error "~S and ~S have different dimensions" ',(first variables) ',var)))
|
|||
|
|
(rest variables))
|
|||
|
|
|
|||
|
|
;; Apply FN with the first two elements (or fewer if size < 2)
|
|||
|
|
(let ((,result (apply ,fn (loop for ,indx below (min ,size 2) collecting
|
|||
|
|
(let ,(map 'list (lambda (var) (list var `(row-major-aref ,var ,indx))) variables)
|
|||
|
|
(progn ,@body))))))
|
|||
|
|
|
|||
|
|
;; Loop over the remaining indices
|
|||
|
|
(loop for ,indx from 2 below ,size do
|
|||
|
|
;; Locally redefine variables to be scalars at a given index
|
|||
|
|
(let ,(mapcar (lambda (var) (list var `(row-major-aref ,var ,indx))) variables)
|
|||
|
|
;; User-supplied function body now evaluated for each index in turn
|
|||
|
|
(setf ,result (funcall ,fn ,result (progn ,@body)))))
|
|||
|
|
,result))))
|
|||
|
|
|
|||
|
|
</code></pre>
|
|||
|
|
<p>[Note: This macro is available in <a href="https://github.com/bendudson/array-operations">this fork</a>
|
|||
|
|
of array-operations, but not Quicklisp]</p>
|
|||
|
|
|
|||
|
|
<p>Using this macro, the maximum value in an array A (of any shape) is:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (vectorize-reduce #'max (a) a)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>The maximum absolute difference between two arrays A and B, of any
|
|||
|
|
shape as long as they have the same shape, is:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (vectorize-reduce #'max (a b) (abs (- a b)))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h2 id="linear-algebra">Linear algebra</h2>
|
|||
|
|
|
|||
|
|
<p>Several packages provide bindings to BLAS and LAPACK libraries,
|
|||
|
|
including:</p>
|
|||
|
|
|
|||
|
|
<ul>
|
|||
|
|
<li><a href="https://github.com/tpapp/lla">lla</a></li>
|
|||
|
|
<li><a href="https://github.com/rigetticomputing/magicl">MAGICL</a></li>
|
|||
|
|
</ul>
|
|||
|
|
|
|||
|
|
<p>A longer list of available packages is on <a href="https://www.cliki.net/linear%20algebra">CLiki’s linear algebra page</a>.</p>
|
|||
|
|
|
|||
|
|
<p>In the examples below the lla package is loaded:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (ql:quickload :lla)
|
|||
|
|
|
|||
|
|
To load "lla":
|
|||
|
|
Load 1 ASDF system:
|
|||
|
|
lla
|
|||
|
|
; Loading "lla"
|
|||
|
|
.
|
|||
|
|
(:LLA)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h3 id="matrix-multiplication">Matrix multiplication</h3>
|
|||
|
|
|
|||
|
|
<p>The <a href="https://github.com/tpapp/lla">lla</a> function <code>mm</code> performs
|
|||
|
|
vector-vector, matrix-vector and matrix-matrix
|
|||
|
|
multiplication.</p>
|
|||
|
|
|
|||
|
|
<h4 id="vector-dot-product">Vector dot product</h4>
|
|||
|
|
|
|||
|
|
<p>Note that one vector is treated as a row vector, and the other as
|
|||
|
|
column:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (lla:mm #(1 2 3) #(2 3 4))
|
|||
|
|
20
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h4 id="matrix-vector-product">Matrix-vector product</h4>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (lla:mm #2A((1 1 1) (2 2 2) (3 3 3)) #(2 3 4))
|
|||
|
|
#(9.0d0 18.0d0 27.0d0)
|
|||
|
|
</code></pre>
|
|||
|
|
<p>which has performed the sum over <code>j</code> of <code>A[i j] * x[j]</code></p>
|
|||
|
|
|
|||
|
|
<h4 id="matrix-matrix-multiply">Matrix-matrix multiply</h4>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (lla:mm #2A((1 2 3) (1 2 3) (1 2 3)) #2A((2 3 4) (2 3 4) (2 3 4)))
|
|||
|
|
#2A((12.0d0 18.0d0 24.0d0) (12.0d0 18.0d0 24.0d0) (12.0d0 18.0d0 24.0d0))
|
|||
|
|
</code></pre>
|
|||
|
|
<p>which summed over <code>j</code> in <code>A[i j] * B[j k]</code></p>
|
|||
|
|
|
|||
|
|
<p>Note that the type of the returned arrays are simple arrays,
|
|||
|
|
specialised to element type <code>double-float</code></p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (type-of (lla:mm #2A((1 0 0) (0 1 0) (0 0 1)) #(1 2 3)))
|
|||
|
|
(SIMPLE-ARRAY DOUBLE-FLOAT (3))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h4 id="outer-product">Outer product</h4>
|
|||
|
|
|
|||
|
|
<p>The <a href="https://github.com/tpapp/array-operations">array-operations</a>
|
|||
|
|
package contains a generalised <a href="https://en.wikipedia.org/wiki/Outer_product">outer product</a>
|
|||
|
|
function:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (ql:quickload :array-operations)
|
|||
|
|
To load "array-operations":
|
|||
|
|
Load 1 ASDF system:
|
|||
|
|
array-operations
|
|||
|
|
; Loading "array-operations"
|
|||
|
|
|
|||
|
|
(:ARRAY-OPERATIONS)
|
|||
|
|
* (aops:outer #'* #(1 2 3) #(2 3 4))
|
|||
|
|
#2A((2 3 4) (4 6 8) (6 9 12))
|
|||
|
|
</code></pre>
|
|||
|
|
<p>which has created a new 2D array <code>A[i j] = B[i] * C[j]</code>. This <code>outer</code>
|
|||
|
|
function can take an arbitrary number of inputs, and inputs with
|
|||
|
|
multiple dimensions.</p>
|
|||
|
|
|
|||
|
|
<h3 id="matrix-inverse">Matrix inverse</h3>
|
|||
|
|
|
|||
|
|
<p>The direct inverse of a dense matrix can be calculated with <code>invert</code></p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (lla:invert #2A((1 0 0) (0 1 0) (0 0 1)))
|
|||
|
|
#2A((1.0d0 0.0d0 -0.0d0) (0.0d0 1.0d0 -0.0d0) (0.0d0 0.0d0 1.0d0))
|
|||
|
|
</code></pre>
|
|||
|
|
<p>e.g</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
|
|||
|
|
A
|
|||
|
|
* (defparameter b (lla:invert a))
|
|||
|
|
B
|
|||
|
|
* (lla:mm a b)
|
|||
|
|
#2A((1.0d0 2.220446049250313d-16 0.0d0)
|
|||
|
|
(0.0d0 1.0d0 0.0d0)
|
|||
|
|
(0.0d0 1.1102230246251565d-16 0.9999999999999998d0))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Calculating the direct inverse is generally not advisable,
|
|||
|
|
particularly for large matrices. Instead the
|
|||
|
|
<a href="https://en.wikipedia.org/wiki/LU_decomposition">LU decomposition</a> can be
|
|||
|
|
calculated and used for multiple inversions.</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
|
|||
|
|
A
|
|||
|
|
* (defparameter b (lla:mm a #(1 2 3)))
|
|||
|
|
B
|
|||
|
|
* (lla:solve (lla:lu a) b)
|
|||
|
|
#(1.0d0 2.0d0 3.0d0)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h3 id="singular-value-decomposition">Singular value decomposition</h3>
|
|||
|
|
|
|||
|
|
<p>The <code>svd</code> function calculates the <a href="https://en.wikipedia.org/wiki/Singular-value_decomposition">singular value decomposition</a>
|
|||
|
|
of a given matrix, returning an object with slots for the three returned
|
|||
|
|
matrices:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
|
|||
|
|
A
|
|||
|
|
* (defparameter a-svd (lla:svd a))
|
|||
|
|
A-SVD
|
|||
|
|
* a-svd
|
|||
|
|
#S(LLA:SVD
|
|||
|
|
:U #2A((-0.6494608633564334d0 0.7205486773948702d0 0.24292013188045855d0)
|
|||
|
|
(-0.3744175632000917d0 -0.5810891192666799d0 0.7225973455785591d0)
|
|||
|
|
(-0.6618248071322363d0 -0.3783451320875919d0 -0.6471807210432038d0))
|
|||
|
|
:D #S(CL-NUM-UTILS.MATRIX:DIAGONAL-MATRIX
|
|||
|
|
:ELEMENTS #(5.593122609997059d0 1.2364443401235103d0
|
|||
|
|
0.43380279311714376d0))
|
|||
|
|
:VT #2A((-0.2344460799312531d0 -0.7211054639318696d0 -0.6519524104506949d0)
|
|||
|
|
(0.2767642134809678d0 -0.6924017945853318d0 0.6663192365460215d0)
|
|||
|
|
(-0.9318994611765425d0 -0.02422116311440764d0 0.3619070730398283d0)))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>The diagonal matrix (singular values) and vectors can be accessed with
|
|||
|
|
functions:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">(lla:svd-u a-svd)
|
|||
|
|
#2A((-0.6494608633564334d0 0.7205486773948702d0 0.24292013188045855d0)
|
|||
|
|
(-0.3744175632000917d0 -0.5810891192666799d0 0.7225973455785591d0)
|
|||
|
|
(-0.6618248071322363d0 -0.3783451320875919d0 -0.6471807210432038d0))
|
|||
|
|
|
|||
|
|
* (lla:svd-d a-svd)
|
|||
|
|
#S(CL-NUM-UTILS.MATRIX:DIAGONAL-MATRIX
|
|||
|
|
:ELEMENTS #(5.593122609997059d0 1.2364443401235103d0 0.43380279311714376d0))
|
|||
|
|
|
|||
|
|
* (lla:svd-vt a-svd)
|
|||
|
|
#2A((-0.2344460799312531d0 -0.7211054639318696d0 -0.6519524104506949d0)
|
|||
|
|
(0.2767642134809678d0 -0.6924017945853318d0 0.6663192365460215d0)
|
|||
|
|
(-0.9318994611765425d0 -0.02422116311440764d0 0.3619070730398283d0))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h2 id="matlisp">Matlisp</h2>
|
|||
|
|
|
|||
|
|
<p>The <a href="https://github.com/bharath1097/matlisp/">Matlisp</a> scientific
|
|||
|
|
computation library provides high performance operations on arrays,
|
|||
|
|
including wrappers around BLAS and LAPACK functions.
|
|||
|
|
It can be loaded using quicklisp:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (ql:quickload :matlisp)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>The nickname for <code>matlisp</code> is <code>m</code>. To avoid typing <code>matlisp:</code> or
|
|||
|
|
<code>m:</code> in front of each symbol, you can define your own package which
|
|||
|
|
uses matlisp
|
|||
|
|
(See the <a href="http://www.gigamonkeys.com/book/programming-in-the-large-packages-and-symbols.html">PCL section on packages</a>):</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defpackage :my-new-code
|
|||
|
|
(:use :common-lisp :matlisp))
|
|||
|
|
#<PACKAGE "MY-NEW-CODE">
|
|||
|
|
|
|||
|
|
* (in-package :my-new-code)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>and to use the <code>#i</code> infix reader (note the same name as for
|
|||
|
|
<code>cmu-infix</code>), run:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (named-readtables:in-readtable :infix-dispatch-table)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h3 id="creating-tensors">Creating tensors</h3>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (matlisp:zeros '(2 2))
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
|
|||
|
|
0.000 0.000
|
|||
|
|
0.000 0.000
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Note that by default matrix storage types are <code>double-float</code>.
|
|||
|
|
To create a complex array using <code>zeros</code>, <code>ones</code> and <code>eye</code>, specify the type:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (matlisp:zeros '(2 2) '((complex double-float)))
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: (COMPLEX DOUBLE-FLOAT)>| #(2 2)
|
|||
|
|
0.000 0.000
|
|||
|
|
0.000 0.000
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>As well as <code>zeros</code> and <code>ones</code> there is <code>eye</code> which creates an identity
|
|||
|
|
matrix:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (matlisp:eye '(3 3) '((complex double-float)))
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: (COMPLEX DOUBLE-FLOAT)>| #(3 3)
|
|||
|
|
1.000 0.000 0.000
|
|||
|
|
0.000 1.000 0.000
|
|||
|
|
0.000 0.000 1.000
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h4 id="ranges">Ranges</h4>
|
|||
|
|
|
|||
|
|
<p>To generate 1D arrays there are the <code>range</code> and <code>linspace</code> functions:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (matlisp:range 1 10)
|
|||
|
|
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(9)
|
|||
|
|
1 2 3 4 5 6 7 8 9
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>The <code>range</code> function rounds down it’s final argument to an integer:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (matlisp:range 1 -3.5)
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: SINGLE-FLOAT>| #(5)
|
|||
|
|
1.000 0.000 -1.000 -2.000 -3.000
|
|||
|
|
>
|
|||
|
|
* (matlisp:range 1 3.3)
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: SINGLE-FLOAT>| #(3)
|
|||
|
|
1.000 2.000 3.000
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Linspace is a bit more general, and the values returned include the
|
|||
|
|
end point.</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (matlisp:linspace 1 10)
|
|||
|
|
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(10)
|
|||
|
|
1 2 3 4 5 6 7 8 9 10
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (matlisp:linspace 0 (* 2 pi) 5)
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(5)
|
|||
|
|
0.000 1.571 3.142 4.712 6.283
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Currently <code>linspace</code> requires real inputs, and doesn’t work with complex numbers.</p>
|
|||
|
|
|
|||
|
|
<h4 id="random-numbers-1">Random numbers</h4>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (matlisp:random-uniform '(2 2))
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
|
|||
|
|
0.7287 0.9480
|
|||
|
|
2.6703E-2 0.1834
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">(matlisp:random-normal '(2 2))
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
|
|||
|
|
0.3536 -1.291
|
|||
|
|
-0.3877 -1.371
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>There are functions for other distributions, including
|
|||
|
|
<code>random-exponential</code>, <code>random-beta</code>, <code>random-gamma</code> and
|
|||
|
|
<code>random-pareto</code>.</p>
|
|||
|
|
|
|||
|
|
<h4 id="reader-macros">Reader macros</h4>
|
|||
|
|
|
|||
|
|
<p>The <code>#d</code> and <code>#e</code> reader macros provide a way to create <code>double-float</code>
|
|||
|
|
and <code>single-float</code> tensors:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* #d[1,2,3]
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(3)
|
|||
|
|
1.000 2.000 3.000
|
|||
|
|
>
|
|||
|
|
|
|||
|
|
* #d[[1,2,3],[4,5,6]]
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
|
|||
|
|
1.000 2.000 3.000
|
|||
|
|
4.000 5.000 6.000
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Note that the comma separators are needed.</p>
|
|||
|
|
|
|||
|
|
<h4 id="tensors-from-arrays">Tensors from arrays</h4>
|
|||
|
|
|
|||
|
|
<p>Common lisp arrays can be converted to Matlisp tensors by copying:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (copy #2A((1 2 3)
|
|||
|
|
(4 5 6))
|
|||
|
|
'#.(tensor 'double-float))
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
|
|||
|
|
1.000 2.000 3.000
|
|||
|
|
4.000 5.000 6.000
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Instances of the <code>tensor</code> class can also be created, specifying the
|
|||
|
|
dimensions. The internal storage of <code>tensor</code> objects is a 1D array
|
|||
|
|
(<code>simple-vector</code>) in a slot <code>store</code>.</p>
|
|||
|
|
|
|||
|
|
<p>For example, to create a <code>double-float</code> type tensor:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">(make-instance (tensor 'double-float)
|
|||
|
|
:dimensions (coerce '(2) '(simple-array index-type (*)))
|
|||
|
|
:store (make-array 2 :element-type 'double-float))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h4 id="arrays-from-tensors">Arrays from tensors</h4>
|
|||
|
|
|
|||
|
|
<p>The array store can be accessed using slots:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter vec (m:range 0 5))
|
|||
|
|
* vec
|
|||
|
|
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(5)
|
|||
|
|
0 1 2 3 4
|
|||
|
|
>
|
|||
|
|
* (slot-value vec 'm:store)
|
|||
|
|
#(0 1 2 3 4)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>Multidimensional tensors are also stored in 1D arrays, and are stored
|
|||
|
|
in column-major order rather than the row-major ordering used for
|
|||
|
|
common lisp arrays. A displaced array will therefore be
|
|||
|
|
transposed.</p>
|
|||
|
|
|
|||
|
|
<p>The contents of a tensor can be copied into an array</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (let ((tens (m:ones '(2 3))))
|
|||
|
|
(m:copy tens 'array))
|
|||
|
|
#2A((1.0d0 1.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>or a list:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (m:copy (m:ones '(2 3)) 'cons)
|
|||
|
|
((1.0d0 1.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h3 id="element-access">Element access</h3>
|
|||
|
|
|
|||
|
|
<p>The <code>ref</code> function is the equivalent of <code>aref</code> for standard CL
|
|||
|
|
arrays, and is also setf-able:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (defparameter a (matlisp:ones '(2 3)))
|
|||
|
|
|
|||
|
|
* (setf (ref a 1 1) 2.0)
|
|||
|
|
2.0d0
|
|||
|
|
* a
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
|
|||
|
|
1.000 1.000 1.000
|
|||
|
|
1.000 2.000 1.000
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<h3 id="element-wise-operations-1">Element-wise operations</h3>
|
|||
|
|
|
|||
|
|
<p>The <code>matlisp-user</code> package, loaded when <code>matlisp</code> is loaded, contains
|
|||
|
|
functions for operating element-wise on tensors.</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (matlisp-user:* 2 (ones '(2 3)))
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
|
|||
|
|
2.000 2.000 2.000
|
|||
|
|
2.000 2.000 2.000
|
|||
|
|
>
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
<p>This includes arithmetic operators ‘+’, ‘-‘, ‘*’, ‘/’ and ‘expt’, but
|
|||
|
|
also <code>sqrt</code>,<code>sin</code>,<code>cos</code>,<code>tan</code>, hyperbolic functions, and their inverses.
|
|||
|
|
The <code>#i</code> reader macro recognises many of these, and uses the
|
|||
|
|
<code>matlisp-user</code> functions:</p>
|
|||
|
|
|
|||
|
|
<pre><code class="language-lisp">* (let ((a (ones '(2 2)))
|
|||
|
|
(b (random-normal '(2 2))))
|
|||
|
|
#i( 2 * a + b ))
|
|||
|
|
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
|
|||
|
|
0.9684 3.250
|
|||
|
|
1.593 1.508
|
|||
|
|
>
|
|||
|
|
|
|||
|
|
* (let ((a (ones '(2 2)))
|
|||
|
|
(b (random-normal '(2 2))))
|
|||
|
|
(macroexpand-1 '#i( 2 * a + b )))
|
|||
|
|
(MATLISP-USER:+ (MATLISP-USER:* 2 A) B)
|
|||
|
|
</code></pre>
|
|||
|
|
|
|||
|
|
|
|||
|
|
<p class="page-source">
|
|||
|
|
Page source: <a href="https://github.com/LispCookbook/cl-cookbook/blob/master/arrays.md">arrays.md</a>
|
|||
|
|
</p>
|
|||
|
|
</div>
|
|||
|
|
|
|||
|
|
<script type="text/javascript">
|
|||
|
|
|
|||
|
|
// Don't write the TOC on the index.
|
|||
|
|
if (window.location.pathname != "/cl-cookbook/") {
|
|||
|
|
$("#toc").toc({
|
|||
|
|
content: "#content", // will ignore the first h1 with the site+page title.
|
|||
|
|
headings: "h1,h2,h3,h4"});
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
$("#two-cols + ul").css({
|
|||
|
|
"column-count": "2",
|
|||
|
|
});
|
|||
|
|
$("#contributors + ul").css({
|
|||
|
|
"column-count": "4",
|
|||
|
|
});
|
|||
|
|
</script>
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
<div>
|
|||
|
|
<footer class="footer">
|
|||
|
|
<hr/>
|
|||
|
|
© 2002–2021 the Common Lisp Cookbook Project
|
|||
|
|
</footer>
|
|||
|
|
|
|||
|
|
</div>
|
|||
|
|
<div id="toc-btn">T<br>O<br>C</div>
|
|||
|
|
</div>
|
|||
|
|
|
|||
|
|
<script text="javascript">
|
|||
|
|
HighlightLisp.highlight_auto({className: null});
|
|||
|
|
</script>
|
|||
|
|
|
|||
|
|
<script type="text/javascript">
|
|||
|
|
function duckSearch() {
|
|||
|
|
var searchField = document.getElementById("searchField");
|
|||
|
|
if (searchField && searchField.value) {
|
|||
|
|
var query = escape("site:lispcookbook.github.io/cl-cookbook/ " + searchField.value);
|
|||
|
|
window.location.href = "https://duckduckgo.com/?kj=b2&kf=-1&ko=1&q=" + query;
|
|||
|
|
// https://duckduckgo.com/params
|
|||
|
|
// kj=b2: blue header in results page
|
|||
|
|
// kf=-1: no favicons
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
</script>
|
|||
|
|
|
|||
|
|
<script async defer data-domain="lispcookbook.github.io/cl-cookbook" src="https://plausible.io/js/plausible.js"></script>
|
|||
|
|
|
|||
|
|
</body>
|
|||
|
|
</html>
|