1137 lines
40 KiB
HTML
1137 lines
40 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="icon" href=
|
||
"assets/cl-logo-blue.png"/>
|
||
<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 style="font-size: 120%" href="https://www.udemy.com/course/common-lisp-programming/?couponCode=LISPY-XMAS2023" title="This course is under a paywall on the Udemy platform. Several videos are freely available so you can judge before diving in. vindarel is (I am) the main contributor to this Cookbook."> Discover our contributor's Lisp course with this Christmas coupon.</a> -->
|
||
<!-- <strong> -->
|
||
<!-- Recently added: 18 videos on MACROS. -->
|
||
<!-- </strong> -->
|
||
<!-- <a style="font-size: 90%" href="https://github.com/vindarel/common-lisp-course-in-videos/">Learn more</a>. -->
|
||
<!-- </p> -->
|
||
|
||
<p class="announce">
|
||
📢 New videos: <a href="https://www.youtube.com/watch?v=h_noB1sI_e8">web dev demo part 1</a>, <a href="https://www.youtube.com/watch?v=xnwc7irnc8k">dynamic page with HTMX</a>, <a href="https://www.youtube.com/watch?v=Zpn86AQRVN8">Weblocks demo</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/Lisp-Stat/array-operations">array-operations</a> maintained by @Symbolics 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 is a fork of <a href="https://github.com/bendudson/array-operations">bendudson/array-operations</a> which is a fork of
|
||
<a href="https://github.com/tpapp/array-operations">tpapp/array-operations</a>, the original
|
||
author.</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://lisp-stat.dev/docs/manuals/array-operations/#generate">Array Operations manual on generate</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))
|
||
;; reverse our symbols list,
|
||
;; since we start from the innermost.
|
||
(syms-rev (reverse syms))
|
||
;; innermost dimension first:
|
||
(dims-rev (loop for i from 0 below rank
|
||
collecting (gensym)))
|
||
;; start with innermost expression
|
||
(result `(progn ,@body)))
|
||
;; 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)))
|
||
;; dimensions reversed so that innermost is last:
|
||
(destructuring-bind ,(reverse dims-rev) ,dims
|
||
,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/Lisp-Stat/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/Lisp-Stat/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–2023 the Common Lisp Cookbook Project
|
||
<div>
|
||
📹 Discover <a style="color: darkgrey; text-decoration: underline", href="https://www.udemy.com/course/common-lisp-programming/?referralCode=2F3D698BBC4326F94358">our contributor's Common Lisp video course on Udemy</a>
|
||
</div>
|
||
</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>
|