Back to index/Generated with Clerk from src/emmy/matrix.cljc@f807468
(ns emmy.matrix
"This namespace contains an implementation of a [[Matrix]] datatype and various
operations for creating and working with [[Matrix]] instances.
[[emmy.matrix]] also extends many Emmy generic operations
to the [[Matrix]] datatype."
(:refer-clojure :exclude [get-in some])
(:require [clojure.core :as core]
[emmy.differential :as d]
[emmy.function :as f]
[emmy.generic :as g]
[emmy.polynomial :as poly]
[emmy.series :as series]
[emmy.structure :as s]
[emmy.util :as u]
[emmy.util.aggregate :as ua]
[emmy.value :as v])
#?(:clj
(:import (clojure.lang Associative AFn IFn Sequential))))
(declare fmap identity-like identity? m:=)
#object[emmy.matrix$m_COLON__EQ_ 0x4ffc0ade "
emmy.matrix$m_COLON__EQ_@4ffc0ade"
]
(derive ::square-matrix ::matrix)
nil
(derive ::column-matrix ::matrix)
nil
(derive ::row-matrix ::matrix)
nil
(derive ::matrix ::f/cofunction)
nil
(deftype Matrix [r c v]
v/IKind
(kind [_] (cond (= r c) ::square-matrix
(= r 1) ::row-matrix
(= c 1) ::column-matrix
:else ::matrix))
d/IPerturbed
(perturbed? [_] (boolean (core/some d/perturbed? v)))
(replace-tag [M old new] (fmap #(d/replace-tag % old new) M))
(extract-tangent [M tag] (fmap #(d/extract-tangent % tag) M))
f/IArity
(arity [_] (transduce (map f/seq-arity) f/combine-arities v))
#?@(:clj
[Object
(equals [this that] (m:= this that))
(toString [_] (pr-str v))
Sequential
Associative
(assoc [_ k entry] (Matrix. r c (assoc v k entry)))
(containsKey [_ k] (contains? v k))
(entryAt [_ k] (.entryAt ^Associative v k))
(count [_] (count v))
(seq [_] (seq v))
(valAt [_ key] (get v key))
(valAt [_ key default] (get v key default))
(empty [this] (fmap g/zero-like this))
(equiv [this that] (m:= this that))
IFn
(invoke [_ a]
(Matrix. r c (mapv (fn [row] (mapv #(% a) row)) v)))
(invoke [_ a b]
(Matrix. r c (mapv (fn [row] (mapv #(% a b) row)) v)))
(invoke [_ a b cx]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx) row)) v)))
(invoke [_ a b cx d]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d) row)) v)))
(invoke [_ a b cx d e]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e) row)) v)))
(invoke [_ a b cx d e f]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f) row)) v)))
(invoke [_ a b cx d e f g]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g) row)) v)))
(invoke [_ a b cx d e f g h]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h) row)) v)))
(invoke [_ a b cx d e f g h i]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i) row)) v)))
(invoke [_ a b cx d e f g h i j]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j) row)) v)))
(invoke [_ a b cx d e f g h i j k]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k) row)) v)))
(invoke [_ a b cx d e f g h i j k l]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l) row)) v)))
(invoke [_ a b cx d e f g h i j k l m]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m) row)) v)))
(invoke [_ a b cx d e f g h i j k l m n]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n) row)) v)))
(invoke [_ a b cx d e f g h i j k l m n o]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o) row)) v)))
(invoke [_ a b cx d e f g h i j k l m n o p]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o p) row)) v)))
(invoke [_ a b cx d e f g h i j k l m n o p q]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o p q) row)) v)))
(invoke [_ a b cx d e f g h i j k l m n o p q rx]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o p q rx) row)) v)))
(invoke [_ a b cx d e f g h i j k l m n o p q rx s]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o p q rx s) row)) v)))
(invoke [_ a b cx d e f g h i j k l m n o p q rx s t]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o p q rx s t) row)) v)))
(invoke [_ a b cx d e f g h i j k l m n o p q rx s t rest]
(Matrix. r c (mapv (fn [row] (mapv #(apply % a b cx d e f g h i j k l m n o p q rx s t rest) row)) v)))
(applyTo [m xs]
(AFn/applyToHelper m xs))]
:cljs
[Object
(toString [_] (pr-str v))
IPrintWithWriter
(-pr-writer [x writer _]
(write-all writer
"#object[emmy.structure.Matrix \""
(.toString x)
"\"]"))
IEmptyableCollection
(-empty [this] (g/zero-like this))
ISequential
IEquiv
(-equiv [this that] (m:= this that))
ISeqable
(-seq [_] (-seq v))
ICounted
(-count [_] (-count v))
IIndexed
(-nth [_ n] (-nth v n))
(-nth [_ n not-found] (-nth v n not-found))
ILookup
(-lookup [_ k] (-lookup v k))
(-lookup [_ k not-found] (-lookup v k not-found))
IAssociative
(-assoc [_ k entry] (Matrix. r c (-assoc v k entry)))
(-contains-key? [_ k] (-contains-key? v k))
IFind
(-find [_ n] (-find v n))
IFn
(-invoke [_ a]
(Matrix. r c (mapv (fn [row] (mapv #(% a) row)) v)))
(-invoke [_ a b]
(Matrix. r c (mapv (fn [row] (mapv #(% a b) row)) v)))
(-invoke [_ a b cx]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx) row)) v)))
(-invoke [_ a b cx d]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d) row)) v)))
(-invoke [_ a b cx d e]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e) row)) v)))
(-invoke [_ a b cx d e f]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f) row)) v)))
(-invoke [_ a b cx d e f g]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g) row)) v)))
(-invoke [_ a b cx d e f g h]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h) row)) v)))
(-invoke [_ a b cx d e f g h i]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i) row)) v)))
(-invoke [_ a b cx d e f g h i j]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j) row)) v)))
(-invoke [_ a b cx d e f g h i j k]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k) row)) v)))
(-invoke [_ a b cx d e f g h i j k l]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l) row)) v)))
(-invoke [_ a b cx d e f g h i j k l m]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m) row)) v)))
(-invoke [_ a b cx d e f g h i j k l m n]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n) row)) v)))
(-invoke [_ a b cx d e f g h i j k l m n o]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o) row)) v)))
(-invoke [_ a b cx d e f g h i j k l m n o p]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o p) row)) v)))
(-invoke [_ a b cx d e f g h i j k l m n o p q]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o p q) row)) v)))
(-invoke [_ a b cx d e f g h i j k l m n o p q rx]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o p q rx) row)) v)))
(-invoke [_ a b cx d e f g h i j k l m n o p q rx s]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o p q rx s) row)) v)))
(-invoke [_ a b cx d e f g h i j k l m n o p q rx s t]
(Matrix. r c (mapv (fn [row] (mapv #(% a b cx d e f g h i j k l m n o p q rx s t) row)) v)))
(-invoke [_ a b cx d e f g h i j k l m n o p q rx s t rest]
(Matrix. r c (mapv (fn [row] (mapv #(apply % a b cx d e f g h i j k l m n o p q rx s t rest) row)) v)))]))
#object[emmy.matrix$eval86072$__GT_Matrix__86162 0x3e4896f7 "
emmy.matrix$eval86072$__GT_Matrix__86162@3e4896f7"
]
(defn matrix?
"Returns true if the supplied `m` is an instance of [[Matrix]], false
otherwise."
[m]
(instance? Matrix m))
#object[emmy.matrix$matrix_QMARK_ 0x7e016978 "
emmy.matrix$matrix_QMARK_@7e016978"
]
(defn num-rows
"Returns the number of rows of the supplied matrix `m`. Throws if a
non-matrix is supplied."
[m]
(if (matrix? m)
(.-r ^Matrix m)
(u/illegal (str "non-matrix supplied: " m))))
#object[emmy.matrix$num_rows 0x4e412c66 "
emmy.matrix$num_rows@4e412c66"
]
(defn num-cols
"Returns the number of columns of the supplied matrix `m`. Throws if a
non-matrix is supplied."
[m]
(if (matrix? m)
(.-c ^Matrix m)
(u/illegal (str "non-matrix supplied: " m))))
#object[emmy.matrix$num_cols 0x61c4658a "
emmy.matrix$num_cols@61c4658a"
]
(defn matrix->vector
"If `m` is already a vector, acts as identity. Else, returns the matrix as a
vector of rows (or throws if neither of these types is passed)."
[m]
(cond (vector? m) m
(matrix? m) (.-v ^Matrix m)
:else
(u/illegal (str "non-matrix supplied: " m))))
#object[emmy.matrix$matrix__GT_vector 0x1bc1a86b "
emmy.matrix$matrix__GT_vector@1bc1a86b"
]
(defn square?
"Returns true if `m` is a square matrix, false otherwise."
[m]
(and (matrix? m)
(= (num-rows m)
(num-cols m))))
#object[emmy.matrix$square_QMARK_ 0x6efe48c2 "
emmy.matrix$square_QMARK_@6efe48c2"
]
(defn column?
"Returns true if `m` is a matrix with a single column (a 'column matrix'),
false otherwise."
[m]
(and (matrix? m)
(= (num-cols m) 1)))
#object[emmy.matrix$column_QMARK_ 0x291b7b67 "
emmy.matrix$column_QMARK_@291b7b67"
]
(defn row?
"Returns true if `m` is a matrix with a single row (a 'row matrix'), false
otherwise."
[m]
(and (matrix? m)
(= (num-rows m) 1)))
#object[emmy.matrix$row_QMARK_ 0x6b31b3fb "
emmy.matrix$row_QMARK_@6b31b3fb"
]
(defn- m:=
"Returns true if the matrices `this` and `that` are of identical shape and
return `v/=` for all entries, false otherwise."
[^Matrix this that]
(and (instance? Matrix that)
(let [^Matrix m that]
(and (= (.-r this) (.-r m))
(= (.-c this) (.-c m))
(v/= (.-v this) (.-v m))))))
#object[emmy.matrix$m_COLON__EQ_ 0x4ffc0ade "
emmy.matrix$m_COLON__EQ_@4ffc0ade"
]
(declare make-diagonal)
#object[emmy.matrix$make_diagonal 0x7da5cc06 "
emmy.matrix$make_diagonal@7da5cc06"
]
(defn- matrix=scalar
"Returns true if the matrix `m` is a diagonal matrix with all diagonal elements
equal to `c`, false otherwise."
[m c]
(and (square? m)
(m:= m (make-diagonal (num-rows m) c))))
#object[emmy.matrix$matrix_EQ_scalar 0x61812328 "
emmy.matrix$matrix_EQ_scalar@61812328"
]
(defn- scalar=matrix
"Returns true if the matrix `m` is a diagonal matrix with all diagonal elements
equal to `c`, false otherwise."
[c m]
(and (square? m)
(m:= (make-diagonal (num-rows m) c) m)))
#object[emmy.matrix$scalar_EQ_matrix 0x1333bc3c "
emmy.matrix$scalar_EQ_matrix@1333bc3c"
]
(defn generate
"Returns a matrix with `r` rows and `c` columns, whose entries are generated by
the supplied function `f`.
If you only supply one dimension `n` the returned matrix will be square.
The entry in the `i`th row and `j`-th column is `(f i j)`."
([n f]
(generate n n f))
([r c f]
(->Matrix r c
(mapv (fn [i]
(mapv (fn [j]
(f i j))
(range c)))
(range r)))))
#object[emmy.matrix$generate 0x757b6bca "
emmy.matrix$generate@757b6bca"
]
(defn literal-matrix
"Generates a `nrows` x `ncols` matrix of symbolic entries, each prefixed by the
supplied symbol `sym`.
If `ncols` (the third argument) is not supplied, returns a square matrix of
size `nrows` x `nrows`.
NOTE: The symbols in the returned matrix record their Einstein-notation path
into the structure that this matrix represents; a `down` of `up` columns. This
means that the returned indices embedded in the symbols look flipped, `ji` vs
`ij`.
For example:
```clojure
(= (literal-matrix 'x 2 2)
(by-rows ['x_0↑0 'x_1↑0]
['x_0↑1 'x_1↑1]))
```"
([sym nrows]
(literal-matrix sym nrows nrows))
([sym nrows ncols]
(let [prefix (str sym "_")]
(generate nrows ncols
(fn [i j]
(symbol (str prefix j "↑" i)))))))
#object[emmy.matrix$literal_matrix 0x26ed5c9f "
emmy.matrix$literal_matrix@26ed5c9f"
]
(defn literal-column-matrix
"Returns a column matrix of `nrows` symbolic entries, each prefixed by the
supplied symbol `sym`.
For example:
```clojure
(= (literal-column-matrix 'x 3)
(by-cols ['x↑0 'x↑1 'x↑2]))
```"
[sym nrows]
(generate nrows 1
(fn [i _]
(symbol
(str sym "↑" i)))))
#object[emmy.matrix$literal_column_matrix 0x4e9da309 "
emmy.matrix$literal_column_matrix@4e9da309"
]
(defn literal-row-matrix
"Returns a row matrix of `ncols` symbolic entries, each prefixed by the
supplied symbol `sym`.
For example:
```clojure
(= (literal-row-matrix 'x 3)
(by-rows ['x_0 'x_1 'x_2]))
```"
[sym ncols]
(generate 1 ncols
(fn [_ j]
(symbol
(str sym "_" j)))))
#object[emmy.matrix$literal_row_matrix 0x404134ba "
emmy.matrix$literal_row_matrix@404134ba"
]
(defn get-in
"Like [[clojure.core/get-in]] for matrices, but obeying the scmutils convention:
only one index is required to get an unboxed element from a column vector.
NOTE that this is perhaps an unprincipled exception..."
[m is]
(let [e (core/get-in m is)]
(if (and (column? m)
(= 1 (count is)))
(e 0)
e)))
#object[emmy.matrix$get_in 0x6eca3243 "
emmy.matrix$get_in@6eca3243"
]
(defn some
"Returns true if `f` is true for some element of the matrix `m`, false
otherwise. (Also works on arbitrary nested sequences.)"
[f m]
(core/some f (flatten m)))
#object[emmy.matrix$some 0x228f92a1 "
emmy.matrix$some@228f92a1"
]
(defn fmap
"Maps `f` over the elements of the matrix `m` returning a new matrix of the same
dimensions as `m`."
[f m]
(->Matrix (num-rows m)
(num-cols m)
(mapv #(mapv f %) m)))
#object[emmy.matrix$fmap 0x129f812 "
emmy.matrix$fmap@129f812"
]
(defn fmap-indexed
"Maps `f` over three arguments:
- each element of the matrix `m`
- its row `i`
- its column `j`
and returns a new matrix of the same dimensions as `m`. "
[f m]
(letfn [(process-row [i row]
(into [] (map-indexed
(fn [j elem] (f elem i j))
row)))]
(let [new-rows (into [] (map-indexed
process-row m))]
(->Matrix (num-rows m)
(num-cols m)
new-rows))))
#object[emmy.matrix$fmap_indexed 0x4804f869 "
emmy.matrix$fmap_indexed@4804f869"
]
(defn- well-formed?
"Returns true if the supplied sequence contains only sequences of the same
length (that could be transformed into columns of a matrix), false otherwise"
[vs]
{:pre [(seq vs) (every? seq vs)]}
(let [counts (map count vs)]
(every? #(= % (first counts))
(next counts))))
#object[emmy.matrix$well_formed_QMARK_ 0x77639300 "
emmy.matrix$well_formed_QMARK_@77639300"
]
(defn by-rows*
"Returns a matrix whose rows consist of the supplied sequence of `rows`. These
all must be the same length.
for a variadic equivalent, see [[by-rows]]."
[rows]
(if (well-formed? rows)
(->Matrix (count rows)
(count (first rows))
(mapv vec rows))
(u/illegal "malformed matrix")))
#object[emmy.matrix$by_rows_STAR_ 0x29270dd7 "
emmy.matrix$by_rows_STAR_@29270dd7"
]
(defn by-rows
"Returns a matrix whose rows consist of the supplied sequence of `rows`. These
all must be the same length.
Variadic equivalent to [[by-rows*]]."
[& rows]
(by-rows* rows))
#object[emmy.matrix$by_rows 0x67809989 "
emmy.matrix$by_rows@67809989"
]
(defn by-cols*
"Returns a matrix whose columns consist of the supplied sequence of `cols`.
These all must be the same length.
for a variadic equivalent, see [[by-cols]]."
[cols]
(if (well-formed? cols)
(->Matrix (count (first cols))
(count cols)
(apply mapv vector cols))
(u/illegal "malformed matrix")))
#object[emmy.matrix$by_cols_STAR_ 0x8a8f2a0 "
emmy.matrix$by_cols_STAR_@8a8f2a0"
]
(defn by-cols
"Returns a matrix whose columns consist of the supplied sequence of `cols`.
These all must be the same length.
Variadic equivalent to [[by-cols*]]."
[& cols]
(by-cols* cols))
#object[emmy.matrix$by_cols 0x63b89b26 "
emmy.matrix$by_cols@63b89b26"
]
(defn row*
"Returns a row matrix populated by the supplied `xs`. For a variadic equivalent,
see [[row]]."
[xs]
{:pre [(not-empty xs)]}
(->Matrix 1 (count xs) [(vec xs)]))
#object[emmy.matrix$row_STAR_ 0x72607b33 "
emmy.matrix$row_STAR_@72607b33"
]
(defn row
"Returns a row matrix populated by the supplied `xs`. Variadic equivalent
to [[row*]]."
[& xs]
(row* xs))
#object[emmy.matrix$row 0x11cfe827 "
emmy.matrix$row@11cfe827"
]
(defn column*
"Returns a column matrix populated by the supplied `xs`. For a variadic equivalent,
see [[column]]."
[xs]
{:pre [(not-empty xs)]}
(->Matrix (count xs) 1 (mapv vector xs)))
#object[emmy.matrix$column_STAR_ 0x5ebd354c "
emmy.matrix$column_STAR_@5ebd354c"
]
(defn column
"Returns a column matrix populated by the supplied `xs`. Variadic equivalent
to [[column*]]."
[& xs]
(column* xs))
#object[emmy.matrix$column 0x79fbe416 "
emmy.matrix$column@79fbe416"
]
(defn transpose
"Returns the transpose of the matrix `m`. The transpose is the original matrix,
with rows and columns swapped."
[m]
(generate (num-cols m)
(num-rows m)
#(core/get-in m [%2 %1])))
#object[emmy.matrix$transpose 0x833a1f9 "
emmy.matrix$transpose@833a1f9"
]
(defn ->structure
"Returns a structure generated by converting `m` into a nested structure with
the supplied `outer-orientation` and `inner-orientation`.
If `t?` is true, the columns of `m` will form the inner tuples. If `t?` is
false, the rows of `m` will form the inner tuples.
By default, if you supply a single argument (the matrix `m`), a matrix turns
into a single outer `::s/down` of inner columns represented as `::up`
structures."
([m] (->structure m ::s/down ::s/up true))
([m outer-orientation inner-orientation t?]
{:pre [(s/valid-orientation? outer-orientation)
(s/valid-orientation? inner-orientation)]}
(let [m' (if t? (transpose m) m)]
(s/make outer-orientation
(mapv #(s/make inner-orientation %)
m')))))
#object[emmy.matrix$__GT_structure 0x6653ee08 "
emmy.matrix$__GT_structure@6653ee08"
]
(defn seq->
"Convert a sequence `xs` (typically, of function arguments) to an up-structure.
Any matrix present in the argument list will be converted to row of columns
via [[->structure]]."
[xs]
(s/up* (map (fn [m]
(if (matrix? m)
(->structure m)
m))
xs)))
#object[emmy.matrix$seq__GT_ 0x5173a0db "
emmy.matrix$seq__GT_@5173a0db"
]
(defn- mul
"Returns the matrix product of `a` and `b`. Throws if `a` and `b` are
incompatible for multiplication."
[a b]
(let [ra (num-rows a)
rb (num-rows b)
ca (num-cols a)
cb (num-cols b)]
(when-not (= ca rb)
(u/illegal "matrices incompatible for multiplication"))
(generate ra cb #(reduce
g/+ (for [k (range ca)]
(g/* (core/get-in a [%1 k])
(core/get-in b [k %2])))))))
#object[emmy.matrix$mul 0x778e25b7 "
emmy.matrix$mul@778e25b7"
]
(defn- scalar*matrix
"Returns a matrix of the same dimensions as matrix `m` with each entry
multiplied (on the left) by the scalar quantity `c`."
[c m]
(fmap #(g/* c %) m))
#object[emmy.matrix$scalar_STAR_matrix 0x62a77ed6 "
emmy.matrix$scalar_STAR_matrix@62a77ed6"
]
(defn- matrix*scalar
"Returns a matrix of the same dimensions as matrix `m` with each entry
multiplied (on the right) by the scalar quantity `c`."
[m c]
(fmap #(g/* % c) m))
#object[emmy.matrix$matrix_STAR_scalar 0x1b7dff24 "
emmy.matrix$matrix_STAR_scalar@1b7dff24"
]
(defn ^:no-doc elementwise
"Applies `f` elementwise between the matrices `a` and `b`. Throws if `a` and `b`
don't have the same dimensions."
[f a b]
(let [ra (num-rows a)
rb (num-rows b)
ca (num-cols a)
cb (num-cols b)]
(when (or (not= ra rb)
(not= ca cb))
(u/illegal "matrices incompatible for operation"))
(generate ra ca #(f (core/get-in a [%1 %2])
(core/get-in b [%1 %2])))))
#object[emmy.matrix$elementwise 0x57c4e1e8 "
emmy.matrix$elementwise@57c4e1e8"
]
(defn two-tensor->
"Converts the square structure `s` into a matrix, and calls the supplied
continuation `cont` with
- the generated matrix
- a function which will restore a matrix to a structure with the same inner
and outer orientations as s
Returns the result of the continuation call."
[s cont]
(if-let [{:keys [inner-size
outer-size
inner-orientation
outer-orientation]}
(s/two-tensor-info s)]
(let [transpose? (= inner-orientation ::s/up)
[major-size minor-size]
(if transpose?
[inner-size outer-size]
[outer-size inner-size])
M (generate major-size minor-size
(fn [i j]
(let [path (if transpose? [j i] [i j])]
(core/get-in s path))))
restore-fn (fn [m]
(->structure
m outer-orientation inner-orientation transpose?))]
(cont M restore-fn))
(u/illegal (str "structure " s " is not a 2-tensor"))))
#object[emmy.matrix$two_tensor__GT_ 0x5161d53a "
emmy.matrix$two_tensor__GT_@5161d53a"
]
(defn two-tensor-operation
"Applies matrix operation `f` to square structure `s` and returns a structure of
the same type as the supplied structure."
[s f]
(two-tensor-> s (fn [m ->s] (->s (f m)))))
#object[emmy.matrix$two_tensor_operation 0x869400e "
emmy.matrix$two_tensor_operation@869400e"
]
(defn structure->matrix
"Given some 2-tensor-shaped structure `s`, returns the corresponding matrix.
The outer orientation is ignored; If the inner structures are `up`, they're
treated as columns. Inner `down` structures are treated as rows."
[s]
(two-tensor-> s (fn [m _] m)))
#object[emmy.matrix$structure__GT_matrix 0x7c991af8 "
emmy.matrix$structure__GT_matrix@7c991af8"
]
(defn- M*u
"Multiply a matrix by an up structure on the right. The return value is up."
[m u]
(when (not= (num-cols m) (count u))
(u/illegal "matrix and tuple incompatible for multiplication"))
(s/up*
(map (fn [i]
(let [row-i (nth m i)]
(ua/generic-sum
(fn [k]
(g/* (nth row-i k)
(nth u k)))
0 (num-cols m))))
(range (num-rows m)))))
#object[emmy.matrix$M_STAR_u 0xc227f22 "
emmy.matrix$M_STAR_u@c227f22"
]
(defn- d*M
"Multiply a matrix `m` by down tuple `d` on the left. The return value has
orientation `down`."
[d m]
(when (not= (count d) (num-rows m))
(u/illegal "matrix and tuple incompatible for multiplication"))
(s/down*
(map (fn [i]
(ua/generic-sum
(fn [k]
(g/* (get d k)
(core/get-in m [k i])))
0 (num-rows m)))
(range (num-cols m)))))
#object[emmy.matrix$d_STAR_M 0x53b5b399 "
emmy.matrix$d_STAR_M@53b5b399"
]
(def ^:dynamic *careful-conversion*
"Set this dynamic variable to `false` to allow [[s->m]] to operate
on structures for which `(* ls ms rs)` does NOT yield a numerical value."
true)
true
(defn s->m
"Convert the structure `ms`, which would be a scalar if the (compatible)
multiplication `(* ls ms rs)` were performed, to a matrix."
([ms rs]
(let [ls (s/compatible-shape (g/* ms rs))]
(s->m ls ms rs)))
([ls ms rs]
(when *careful-conversion*
(assert (v/numerical? (g/* ls (g/* ms rs)))))
(let [ndowns (s/dimension ls)
nups (s/dimension rs)]
(generate ndowns nups
(fn [i j]
(g/* (s/unflatten
(s/basis-unit ndowns i) ls)
(g/* ms (s/unflatten
(s/basis-unit nups j) rs))))))))
#object[emmy.matrix$s__GT_m 0x119ee0d1 "
emmy.matrix$s__GT_m@119ee0d1"
]
(defn as-matrix
"Any one argument function of a structure can be seen as a matrix. This is only
useful if the function has a linear multiplier (e.g. derivative)"
[F]
(fn [s]
(let [v (F s)]
(s->m v s))))
#object[emmy.matrix$as_matrix 0x256476c0 "
emmy.matrix$as_matrix@256476c0"
]
(defn nth-row
"Returns the `n`-th row of the supplied matrix `m` as a `down` structure."
[m n]
(s/down* (get m n)))
#object[emmy.matrix$nth_row 0x3aa0d621 "
emmy.matrix$nth_row@3aa0d621"
]
(defn nth-col
"Returns the `n`-th column of the supplied matrix `m` as an `up` structure."
[m n]
(s/up* (map #(% n) m)))
#object[emmy.matrix$nth_col 0x44779397 "
emmy.matrix$nth_col@44779397"
]
(defn diagonal
"Returns the diagonal of the supplied matrix `m` as an up structure. Errors if a
type other than a diagonal matrix is supplied."
[m]
{:pre [(square? m)]}
(let [rows (num-rows m)]
(s/up* (map #(core/get-in m [% %])
(range 0 rows)))))
#object[emmy.matrix$diagonal 0x22096ba0 "
emmy.matrix$diagonal@22096ba0"
]
(defn up->column-matrix
"Returns a column matrix with the contents of the supplied `up` structure.
Errors if any other type is provided."
[v]
{:pre [(s/up? v)]}
(column* v))
#object[emmy.matrix$up__GT_column_matrix 0x245ddb57 "
emmy.matrix$up__GT_column_matrix@245ddb57"
]
(defn column-matrix->up
"Returns the single column from the supplied column matrix as an `up`. Errors if
some other type is supplied."
[m]
{:pre [(column? m)]}
(nth-col m 0))
#object[emmy.matrix$column_matrix__GT_up 0x18f01023 "
emmy.matrix$column_matrix__GT_up@18f01023"
]
(defn column-matrix->vector
"Returns the single column from the supplied column matrix as a vector. Errors
if some other type is supplied."
[m]
{:pre [(column? m)]}
(mapv first m))
#object[emmy.matrix$column_matrix__GT_vector 0x60e3c381 "
emmy.matrix$column_matrix__GT_vector@60e3c381"
]
(defn down->row-matrix
"Returns a row matrix with the contents of the supplied `down` structure.
Errors if any other type is provided."
[v]
{:pre [(s/down? v)]}
(by-rows (s/structure->vector v)))
#object[emmy.matrix$down__GT_row_matrix 0x95ef056 "
emmy.matrix$down__GT_row_matrix@95ef056"
]
(defn row-matrix->down
"Returns the single row from the supplied row matrix as a `down`. Errors if some
other type is supplied."
[m]
{:pre [(row? m)]}
(nth-row m 0))
#object[emmy.matrix$row_matrix__GT_down 0x43130549 "
emmy.matrix$row_matrix__GT_down@43130549"
]
(defn row-matrix->vector
"Returns the single row from the supplied row matrix as a vector. Errors if some
other type is supplied."
[m]
{:pre [(row? m)]}
(nth m 0))
#object[emmy.matrix$row_matrix__GT_vector 0x5b7a7400 "
emmy.matrix$row_matrix__GT_vector@5b7a7400"
]
(defn m->s
"Convert the matrix `m` into a structure `S`, guided by the requirement that `(*
ls S rs)` should be a scalar."
[ls m rs]
(let [ncols (num-cols m)
col-shape (s/compatible-shape ls)
ms (s/unflatten (for [j (range ncols)]
(s/unflatten (nth-col m j) col-shape))
(s/compatible-shape rs))]
(when *careful-conversion*
(assert (v/numerical? (g/* ls (g/* ms rs)))
(str "product is not numerical: " ls ms rs)))
ms))
#object[emmy.matrix$m__GT_s 0x25fb4b25 "
emmy.matrix$m__GT_s@25fb4b25"
]
(defn s:transpose
"Given structural inputs `ls` (optional), `ms` and `rs`, constrained such
that `(* ls ms rs)` returns a numerical quantity, returns a result such that
the following relationship remains true:
```clj
(* <ls| (* ms |rs>)) = (* <rs| (* (s:transpose ms) |ls>))
```
For example:
```clj
(let [ls (s/up 1 2)
ms (s/up (s/down 1 2) (s/down 3 4))
rs (s/down 1 2)]
(g/* ls (g/* ms rs))
;;=> 27
(g/* rs (g/* (s:transpose ls ms rs) ls))
;;=> 27
)
```
`ls` is optional. If `ls` is not supplied, a compatible shape is generated
internally."
([ms rs]
(let [ls (s/compatible-shape (g/* ms rs))]
(s:transpose ls ms rs)))
([ls ms rs]
(m->s rs (transpose (s->m ls ms rs)) ls)))
#object[emmy.matrix$s_COLON_transpose 0x5b59d38c "
emmy.matrix$s_COLON_transpose@5b59d38c"
]
(defn s:transpose-orientation
"Given some 2 tensor `s`, returns a structure with elements 'transposed' by
swapping the inner and outer orientations and dimensions, like a matrix
transpose.
Orientations are only flipped if they are different in the input. If the inner
and outer orientations of `s` are the same, the returned structure has this
identical orientation.
For example:
```clj
;; opposite orientation gets flipped:
(s:transpose-orientation (s/up (s/down 1 2 3) (s/down 4 5 6)))
;;=> (down (up 1 4) (up 2 5) (up 3 6))
;; same orientation stays the same:
(s:transpose-orientation (s/down (s/down 1 2 3) (s/down 4 5 6)))
;;=> (down (down 1 4) (down 2 5) (down 3 6))
```
See [[structure/two-tensor?]] for more detail on 2 tensors.
NOTE: In scmutils, this function is called `s:transpose2`."
[s]
(let [ret (two-tensor-operation s transpose)]
(if (= (s/orientation ret)
(s/orientation (first ret)))
ret
(s/transpose ret))))
#object[emmy.matrix$s_COLON_transpose_orientation 0xbbace3 "
emmy.matrix$s_COLON_transpose_orientation@bbace3"
]
(declare invert)
#object[emmy.matrix$classical_adjoint_formula$inv__86359 0x8e91280 "
emmy.matrix$classical_adjoint_formula$inv__86359@8e91280"
]
(defn s:invert
"Given some 2-tensor `s` (a 'square' nested structure), returns a structure
that represents the multiplicative inverse of the supplied structure. The
inner and outer structure orientations of `(s:invert s)` are the SAME as `s`.
If `s` is an up-of-downs or down-of-ups, `(g/* s (s:invert s))`
and `(g/* (s:invert s) s)` will evaluate to an identity-matrix-shaped
up-of-downs or down-of-ups.
If `s` is an up-of-ups or down-of-downs, multiplying `s` `(s:invert s)` will
result in a scalar, as both structures collapse.
NOTE: I DO NOT yet understand the meaning of this scalar! If you do, please
open a pull request and explain it here."
[s]
(let [ret (two-tensor-operation s invert)]
(if (= (s/orientation ret)
(s/orientation (first ret)))
(s/transpose ret)
ret)))
#object[emmy.matrix$s_COLON_invert 0x77b4d2db "
emmy.matrix$s_COLON_invert@77b4d2db"
]
(defn- delete
"Returns the vector formed by deleting the `i`'th element of the given vector
`v`."
[v i]
(if (vector? v)
(into (subvec v 0 i)
(subvec v (inc i)))
(delete (into [] v) i)))
#object[emmy.matrix$delete 0x62e4ba12 "
emmy.matrix$delete@62e4ba12"
]
(defn with-substituted-row
"Returns a new matrix of identical shape to `m`, with the vector `v` substituted
for the `i`th row."
[m i v]
{:pre [(matrix? m)
(<= 0 i)
(< i (num-rows m))
(= (num-cols m)
(count v))]}
(assoc m i v))
#object[emmy.matrix$with_substituted_row 0x3f65de19 "
emmy.matrix$with_substituted_row@3f65de19"
]
(defn submatrix
"Returns the submatrix of the matrix (or matrix-like structure) `s` generated by
taking
- rows from `lowrow` -> `hirow` (inclusive)
- columns from `lowcol` -> `hicol` (inclusive)"
[x lowrow hirow lowcol hicol]
(let [m (if (s/structure? x)
(two-tensor-> x (fn [m _] m))
x)]
(generate (inc (- hirow lowrow))
(inc (- hicol lowcol))
(fn [i j]
(core/get-in m [(+ i lowrow)
(+ j lowcol)])))))
#object[emmy.matrix$submatrix 0x7f58b667 "
emmy.matrix$submatrix@7f58b667"
]
(defn without
"Returns the matrix formed by deleting the `i`-th row and `j`-th column of the
given matrix `m`.
This is also called the 'minor' of m."
[m i j]
(->Matrix
(dec (num-rows m))
(dec (num-cols m))
(mapv #(delete % j)
(delete (matrix->vector m) i))))
#object[emmy.matrix$without 0x2e2b915d "
emmy.matrix$without@2e2b915d"
]
(defn- checkerboard-negate [s i j]
(if (even? (+ i j))
s
(g/negate s)))
#object[emmy.matrix$checkerboard_negate 0x69588aa0 "
emmy.matrix$checkerboard_negate@69588aa0"
]
(defn dimension
"Returns the 'dimension', i.e., the number of rows & columns, of the supplied
square matrix. Errors if some other type is supplied."
[m]
{:pre [(square? m)]}
(num-rows m))
#object[emmy.matrix$dimension 0x78dc0484 "
emmy.matrix$dimension@78dc0484"
]
(defn trace
"Returns the trace (the sum of diagonal elements) of the square matrix `m`.
Generic operations are used, so this works on symbolic square matrices."
[m]
{:pre [(square? m)]}
(let [rows (num-rows m)]
(transduce (map #(core/get-in m [% %]))
g/+
(range 0 rows))))
#object[emmy.matrix$trace 0x428d6e76 "
emmy.matrix$trace@428d6e76"
]

Note from GJS in scmutils: "Kleanthes Koniaris determinant routine, slightly edited by GJS"

(defn general-determinant
"Given coefficient procedures `add`, `sub`, `mul` and `zero?`, returns a
procedure that efficiently computes the determinant of the supplied square
matrix `m`.
[[general-determinant]] is useful for generating fast type-specific
determinant routines. See [[determinant]] for a default using generic
arithmetic."
[add sub mul zero?]
(let [zero (add)]
(fn [m]
{:pre [(square? m)]}
(let [c-det (atom nil)]
(letfn [(c-det* [row [col & cols :as active-cols]]
(if-not cols
;; one active column
(core/get-in m [row col])
(loop [idx 0
remaining-cols active-cols
answer zero]
(if-not (seq remaining-cols)
answer
(let [term (core/get-in m [row (first remaining-cols)])]
(if (zero? term)
(recur (inc idx)
(rest remaining-cols)
answer)
(let [without-i (delete active-cols idx)
delta (mul term (@c-det (inc row) without-i))]
(recur (inc idx)
(rest remaining-cols)
(if (even? idx)
(add answer delta)
(sub answer delta))))))))))]
(reset! c-det (memoize c-det*))
(@c-det 0 (range (dimension m))))))))
#object[emmy.matrix$general_determinant 0x305eae03 "
emmy.matrix$general_determinant@305eae03"
]
(def ^{:arglists '([m])}
determinant
"Returns the determinant of the supplied square matrix `m`.
Generic operations are used, so this works on symbolic square matrices."
(general-determinant g/+ g/- g/* g/numeric-zero?))
#object[emmy.matrix$general_determinant$fn__86338 0x1dfe3f22 "
emmy.matrix$general_determinant$fn__86338@1dfe3f22"
]
(defn cofactors
"Returns the matrix of cofactors of the supplied square matrix `m`."
[m]
{:pre [(square? m)]}
(let [r (num-rows m)]
(cond (< r 2) m
(= r 2) (let [[[a b] [c d]] m]
(->Matrix 2 2 [[d (g/negate c)]
[(g/negate b) a]]))
:else (generate r r
(fn [i j]
(-> (without m i j)
(determinant)
(checkerboard-negate i j)))))))
#object[emmy.matrix$cofactors 0x59ee03d "
emmy.matrix$cofactors@59ee03d"
]

The following implements the classical adjoint formula for the inverse of a matrix. This may be useful for symbolic applications.

(defn classical-adjoint-formula
"Given coefficient procedures `add`, `sub`, `mul` and `zero?`, returns a
procedure that efficiently computes the inverse of the supplied square
matrix `m`.
[[classical-adjoint-formula]] is useful for generating fast type-specific
matrix inversion routines. See [[invert]] for a default using generic
arithmetic."
[add sub mul div zero?]
(let [det (general-determinant add sub mul zero?)]
(fn inv [A]
(let [dim (dimension A)]
(if (= dim 1)
(->Matrix 1 1 [[(div (core/get-in A [0 0]))]])
(let* [d (det A)
-d (sub d)]
(generate dim dim
(fn [i j]
(let [denom (if (even? (+ i j)) d -d)]
(div (det (without A j i)) denom))))))))))
#object[emmy.matrix$classical_adjoint_formula 0x2cc9d026 "
emmy.matrix$classical_adjoint_formula@2cc9d026"
]
(def ^{:arglists '([A])}
invert
"Returns the inverse of the supplied square matrix `m`."
(classical-adjoint-formula g/+ g/- g/* g// g/numeric-zero?))
#object[emmy.matrix$classical_adjoint_formula$inv__86359 0x8e91280 "
emmy.matrix$classical_adjoint_formula$inv__86359@8e91280"
]
(defn- m-div-m [m1 m2]
(mul m1 (invert m2)))
#object[emmy.matrix$m_div_m 0x69a24d0a "
emmy.matrix$m_div_m@69a24d0a"
]
(defn- m-div-c
"Returns the result of multiplying (on the right) the scalar `c` by the inverse
of matrix `m`."
[m c]
(matrix*scalar m (g/invert c)))
#object[emmy.matrix$m_div_c 0x5496833 "
emmy.matrix$m_div_c@5496833"
]
(defn- c-div-m
"Returns the result of multiplying (on the left) the scalar `c` by the inverse
of matrix `m`."
[c m]
(scalar*matrix c (invert m)))
#object[emmy.matrix$c_div_m 0x77a4efaf "
emmy.matrix$c_div_m@77a4efaf"
]
(defn s:inverse
([ms rs]
(let [ls (s/compatible-shape (g/* ms rs))]
(s:inverse ls ms rs)))
([ls ms rs]
(m->s (s/compatible-shape rs)
(invert (s->m ls ms rs))
(s/compatible-shape ls))))
#object[emmy.matrix$s_COLON_inverse 0x158e599d "
emmy.matrix$s_COLON_inverse@158e599d"
]
(defn s:solve-linear-left [M product]
(let [cp (s/compatible-shape product)
cr (s/compatible-shape (s/s:* cp M))]
(s/s:* (s:inverse cp M cr) product)))
#object[emmy.matrix$s_COLON_solve_linear_left 0x6e3f9c3b "
emmy.matrix$s_COLON_solve_linear_left@6e3f9c3b"
]
(defn s:solve-linear-right [product M]
(let [cp (s/compatible-shape product)
cr (s/compatible-shape (s/s:* M cp))]
(s/s:* product (s:inverse cr M cp))))
#object[emmy.matrix$s_COLON_solve_linear_right 0x20f9793f "
emmy.matrix$s_COLON_solve_linear_right@20f9793f"
]
(defn s:divide-by-structure [rv s]
(s:solve-linear-left s rv))
#object[emmy.matrix$s_COLON_divide_by_structure 0x70f03318 "
emmy.matrix$s_COLON_divide_by_structure@70f03318"
]
(defn make-zero
"Return a zero-valued matrix of `m` rows and `n` columns (`nXn` if only `n` is
supplied)."
([n] (make-zero n n))
([m n] (generate m n (constantly 0))))
#object[emmy.matrix$make_zero 0x1755de2a "
emmy.matrix$make_zero@1755de2a"
]
(defn I
"Return the identity matrix of order `n`."
[n]
(generate n n s/kronecker))
#object[emmy.matrix$I 0x59946bf "
emmy.matrix$I@59946bf"
]
(defn identity-like
"Return an identity matrix whose ones and zeros match the types of the supplied
square matrix `M`. Errors if a non-square matrix `M` is supplied."
[M]
(if-not (square? M)
(u/illegal "identity-like on non-square")
(fmap-indexed (fn [elem i j]
(if (= i j)
(g/one-like elem)
(g/zero-like elem)))
M)))
#object[emmy.matrix$identity_like 0x6e42f7b2 "
emmy.matrix$identity_like@6e42f7b2"
]
(defn identity?
"Returns true if the supplied matrix `m` is an identity matrix, false
otherwise."
[m]
(and (square? m)
(let [n (dimension m)]
(every? true?
(for [i (range n)
j (range n)
:let [entry (core/get-in m [i j])]]
(if (= i j)
(g/one? entry)
(g/zero? entry)))))))
#object[emmy.matrix$identity_QMARK_ 0x1f5b7191 "
emmy.matrix$identity_QMARK_@1f5b7191"
]
(defn make-diagonal
"Given a single (sequential) argument `v`, returns the diagonal matrix of
order `(count v)` with the elements of the sequence `v` along the diagonal.
Given two arguments `n` and some constant `x`, returns a diagonal `n` by `n`
matrix with `x` in every entry of the diagonal.
`(make-diagonal <n> 1)` is equivalent to `(I n)`."
([v]
(let [v (vec v)
n (count v)]
(generate n n (fn [i j]
(if (= i j) (v i) 0)))))
([n x]
(generate n n #(if (= %1 %2) x 0))))
#object[emmy.matrix$make_diagonal 0x7da5cc06 "
emmy.matrix$make_diagonal@7da5cc06"
]
(defn diagonal?
"Returns true if `m` is a diagonal matrix (i.e., a square matrix where every
non-diagonal element is zero), false otherwise."
[m]
(and (square? m)
(let [n (dimension m)]
(every? true?
(for [i (range n)
j (range n)
:when (not= i j)
:let [entry (core/get-in m [i j])]]
(g/zero? entry))))))
#object[emmy.matrix$diagonal_QMARK_ 0x264648de "
emmy.matrix$diagonal_QMARK_@264648de"
]
(defn symmetric?
"Returns true if the supplied matrix `M` is equal to its own transpose (i.e.,
symmetric), false otherwise."
[M]
(g/zero?
(g/simplify
(g/sub (transpose M) M))))
#object[emmy.matrix$symmetric_QMARK_ 0x69255e17 "
emmy.matrix$symmetric_QMARK_@69255e17"
]
(defn antisymmetric?
"Returns true if the supplied matrix `M` is equal to the negation of its own
transpose (i.e., antisymmetric), false otherwise."
[M]
(g/zero?
(g/simplify
(g/add (transpose M) M))))
#object[emmy.matrix$antisymmetric_QMARK_ 0x1023dbe6 "
emmy.matrix$antisymmetric_QMARK_@1023dbe6"
]
(defn characteristic-polynomial
"Returns the [characteristic
polynomial](https://en.wikipedia.org/wiki/Characteristic_polynomial) of the
square matrix `m`.
If only `m` is supplied, returns a [[polynomial/Polynomial]] instance
representing the matrix `m`'s characteristic polynomial.
If `x` is supplied, returns the value of the characteristic polynomial of `m`
evaluated at `x`.
Typically `x` will be a symbolic variable, but if you wanted to get the value
of the characteristic polynomial at some particular numerical point `x` you
could pass that too."
([m]
(characteristic-polynomial m (poly/identity)))
([m x]
(let [r (num-rows m)
c (num-cols m)]
(when-not (= r c) (u/illegal "not square"))
(let [Ix (make-diagonal r x)]
(determinant
(g/- Ix m))))))
#object[emmy.matrix$characteristic_polynomial 0x54b0f78a "
emmy.matrix$characteristic_polynomial@54b0f78a"
]

Solving

(defn cramers-rule
"Given coefficient procedures `add`, `sub`, `mul`, `div` and `zero?`, returns a
procedure that efficiently computes the solution to an inhomogeneous system of
linear equations, `A*x=b`, where the matrix `A` and the column matrix `b` are
given. The returned procedure returns the column matrix `x`.
Unlike LU decomposition, Cramer's rule generalizes to symbolic solutions.
[[cramers-rule]] is useful for generating fast type-specific linear equation
solvers. See [[solve]] for a default using generic arithmetic."
[add sub mul div zero?]
(let [det (general-determinant add sub mul zero?)]
(fn solve [A b]
{:pre [(square? A)
(column? b)
(= (dimension A)
(num-rows b))]}
(let [bv (nth-col b 0)
bn (num-rows b)
d (det A)
At (transpose A)]
(column*
(mapv (fn [i]
(div (det (with-substituted-row At i bv))
d))
(range bn)))))))
#object[emmy.matrix$cramers_rule 0x40a4457a "
emmy.matrix$cramers_rule@40a4457a"
]
(def ^{:arglists '([A b])}
solve
"Given a matrix `A` and a column matrix `b`, computes the solution
to an inhomogeneous system of linear equations, `A*x=b`, where the matrix `A`
and the column matrix `b` are given.
Returns the column matrix `x`.
Unlike LU decomposition, Cramer's rule generalizes to symbolic solutions."
(cramers-rule g/+ g/- g/* g// g/numeric-zero?))
#object[emmy.matrix$cramers_rule$solve__86436 0x3b27e9f "
emmy.matrix$cramers_rule$solve__86436@3b27e9f"
]
(defn rsolve
"Generalization of [[solve]] that can handle `up` and `down` structures, as well
as `row` and `column` matrices.
Given `row` or `down` values for `b`, `A` is appropriately transposed before
solving."
[b A]
(cond (s/up? b) (column-matrix->up
(solve A (up->column-matrix b)))
(column? b) (solve A b)
(s/down? b) (row-matrix->down
(transpose
(solve (transpose A)
(transpose
(down->row-matrix b)))))
(row? b)
(transpose
(solve (transpose A)
(transpose b)))
:else (u/illegal (str "I don't know how to solve:" b A))))
#object[emmy.matrix$rsolve 0xd4079e4 "
emmy.matrix$rsolve@d4079e4"
]

Generic Operation Installation

(defmethod g/zero? [::matrix] [a] (every? #(every? g/zero? %) a))
#object[clojure.lang.MultiFn 0x78e2923f "
clojure.lang.MultiFn@78e2923f"
]
(defmethod g/one? [::matrix] [_] false)
#object[clojure.lang.MultiFn 0x710fa3bd "
clojure.lang.MultiFn@710fa3bd"
]
(defmethod g/identity? [::matrix] [m] (identity? m))
#object[clojure.lang.MultiFn 0x3e853a24 "
clojure.lang.MultiFn@3e853a24"
]
(defmethod g/zero-like [::matrix] [m] (fmap g/zero-like m))
#object[clojure.lang.MultiFn 0x79337717 "
clojure.lang.MultiFn@79337717"
]
(defmethod g/one-like [::matrix] [m] (identity-like m))
#object[clojure.lang.MultiFn 0x46a24adf "
clojure.lang.MultiFn@46a24adf"
]
(defmethod g/identity-like [::matrix] [m] (identity-like m))
#object[clojure.lang.MultiFn 0x787144b0 "
clojure.lang.MultiFn@787144b0"
]
(defmethod g/freeze [::matrix] [^Matrix m]
(if (= (.-c m) 1)
`(~'column-matrix ~@(map (comp g/freeze first) (.-v m)))
`(~'matrix-by-rows ~@(map #(mapv g/freeze %) (.-v m)))))
#object[clojure.lang.MultiFn 0x7a0db79a "
clojure.lang.MultiFn@7a0db79a"
]
(defmethod g/exact? [::matrix] [^Matrix m] (every? #(every? g/exact? %) (.-v m)))
#object[clojure.lang.MultiFn 0x1864e04e "
clojure.lang.MultiFn@1864e04e"
]
(defmethod v/= [::matrix ::matrix] [a b] (m:= a b))
#object[clojure.lang.MultiFn 0x744f207b "
clojure.lang.MultiFn@744f207b"
]
(defmethod v/= [::square-matrix ::v/scalar] [m c] (matrix=scalar m c))
#object[clojure.lang.MultiFn 0x744f207b "
clojure.lang.MultiFn@744f207b"
]
(defmethod v/= [::v/scalar ::square-matrix] [c m] (scalar=matrix c m))
#object[clojure.lang.MultiFn 0x744f207b "
clojure.lang.MultiFn@744f207b"
]
(defmethod g/negate [::matrix] [a] (fmap g/negate a))
#object[clojure.lang.MultiFn 0x1b9c4031 "
clojure.lang.MultiFn@1b9c4031"
]
(defmethod g/sub [::matrix ::matrix] [a b] (elementwise g/- a b))
#object[clojure.lang.MultiFn 0x3b212ce4 "
clojure.lang.MultiFn@3b212ce4"
]
(defmethod g/sub [::square-matrix ::v/scalar] [a b]
(elementwise g/- a (make-diagonal (num-rows a) b)))
#object[clojure.lang.MultiFn 0x3b212ce4 "
clojure.lang.MultiFn@3b212ce4"
]
(defmethod g/sub [::v/scalar ::square-matrix] [a b]
(elementwise g/- (make-diagonal (num-rows b) a) b))
#object[clojure.lang.MultiFn 0x3b212ce4 "
clojure.lang.MultiFn@3b212ce4"
]
(defmethod g/add [::matrix ::matrix] [a b] (elementwise g/+ a b))
#object[clojure.lang.MultiFn 0x7e1fb1be "
clojure.lang.MultiFn@7e1fb1be"
]
(defmethod g/add [::square-matrix ::v/scalar] [a b]
(elementwise g/+ a (make-diagonal (num-rows a) b)))
#object[clojure.lang.MultiFn 0x7e1fb1be "
clojure.lang.MultiFn@7e1fb1be"
]
(defmethod g/add [::v/scalar ::square-matrix] [a b]
(elementwise g/+ (make-diagonal (num-rows b) a) b))
#object[clojure.lang.MultiFn 0x7e1fb1be "
clojure.lang.MultiFn@7e1fb1be"
]
(defmethod g/mul [::matrix ::matrix] [a b] (mul a b))
#object[clojure.lang.MultiFn 0x4485ec13 "
clojure.lang.MultiFn@4485ec13"
]
(defmethod g/mul [::v/scalar ::matrix] [n a] (scalar*matrix n a))
#object[clojure.lang.MultiFn 0x4485ec13 "
clojure.lang.MultiFn@4485ec13"
]
(defmethod g/mul [::matrix ::v/scalar] [a n] (matrix*scalar a n))
#object[clojure.lang.MultiFn 0x4485ec13 "
clojure.lang.MultiFn@4485ec13"
]
(defmethod g/mul [::matrix ::s/up] [m u] (M*u m u))
#object[clojure.lang.MultiFn 0x4485ec13 "
clojure.lang.MultiFn@4485ec13"
]
(defmethod g/mul [::s/down ::matrix] [d m] (d*M d m))
#object[clojure.lang.MultiFn 0x4485ec13 "
clojure.lang.MultiFn@4485ec13"
]
(defmethod g/div [::matrix ::v/scalar] [m c] (m-div-c m c))
#object[clojure.lang.MultiFn 0x233dab3a "
clojure.lang.MultiFn@233dab3a"
]
(defmethod g/div [::v/scalar ::square-matrix] [c m] (c-div-m c m))
#object[clojure.lang.MultiFn 0x233dab3a "
clojure.lang.MultiFn@233dab3a"
]
(defmethod g/div [::column-matrix ::square-matrix] [c m] (rsolve c m))
#object[clojure.lang.MultiFn 0x233dab3a "
clojure.lang.MultiFn@233dab3a"
]
(defmethod g/div [::row-matrix ::square-matrix] [r m] (rsolve r m))
#object[clojure.lang.MultiFn 0x233dab3a "
clojure.lang.MultiFn@233dab3a"
]
(defmethod g/div [::s/up ::square-matrix] [u m] (rsolve u m))
#object[clojure.lang.MultiFn 0x233dab3a "
clojure.lang.MultiFn@233dab3a"
]
(defmethod g/div [::s/down ::square-matrix] [d m] (rsolve d m))
#object[clojure.lang.MultiFn 0x233dab3a "
clojure.lang.MultiFn@233dab3a"
]
(defmethod g/div [::matrix ::square-matrix] [d m] (m-div-m d m))
#object[clojure.lang.MultiFn 0x233dab3a "
clojure.lang.MultiFn@233dab3a"
]
(defmethod g/exp [::square-matrix] [m] (series/exp-series m))
#object[clojure.lang.MultiFn 0x7168c357 "
clojure.lang.MultiFn@7168c357"
]
(defmethod g/cos [::square-matrix] [m] (series/cos-series m))
#object[clojure.lang.MultiFn 0x3be6105f "
clojure.lang.MultiFn@3be6105f"
]
(defmethod g/sin [::square-matrix] [m] (series/sin-series m))
#object[clojure.lang.MultiFn 0x163caf8e "
clojure.lang.MultiFn@163caf8e"
]
(defmethod g/tan [::square-matrix] [m] (series/tan-series m))
#object[clojure.lang.MultiFn 0x428252da "
clojure.lang.MultiFn@428252da"
]
(defmethod g/sec [::square-matrix] [m] (series/sec-series m))
#object[clojure.lang.MultiFn 0x3b1846b4 "
clojure.lang.MultiFn@3b1846b4"
]
(defmethod g/acos [::square-matrix] [m] (series/acos-series m))
#object[clojure.lang.MultiFn 0x1380d514 "
clojure.lang.MultiFn@1380d514"
]
(defmethod g/asin [::square-matrix] [m] (series/asin-series m))
#object[clojure.lang.MultiFn 0x43c45b34 "
clojure.lang.MultiFn@43c45b34"
]
(defmethod g/atan [::square-matrix] [m] (series/atan-series m))
#object[clojure.lang.MultiFn 0x6a540c36 "
clojure.lang.MultiFn@6a540c36"
]
(defmethod g/acot [::square-matrix] [m] (series/acot-series m))
#object[clojure.lang.MultiFn 0x7c9d3014 "
clojure.lang.MultiFn@7c9d3014"
]
(defmethod g/cosh [::square-matrix] [m] (series/cosh-series m))
#object[clojure.lang.MultiFn 0x8dc9785 "
clojure.lang.MultiFn@8dc9785"
]
(defmethod g/sinh [::square-matrix] [m] (series/sinh-series m))
#object[clojure.lang.MultiFn 0x1d61d948 "
clojure.lang.MultiFn@1d61d948"
]
(defmethod g/tanh [::square-matrix] [m] (series/tanh-series m))
#object[clojure.lang.MultiFn 0x42231f14 "
clojure.lang.MultiFn@42231f14"
]
(defmethod g/asinh [::square-matrix] [m] (series/asinh-series m))
#object[clojure.lang.MultiFn 0x2d117535 "
clojure.lang.MultiFn@2d117535"
]
(defmethod g/atanh [::square-matrix] [m] (series/atanh-series m))
#object[clojure.lang.MultiFn 0x38f36f68 "
clojure.lang.MultiFn@38f36f68"
]
(defmethod g/simplify [::matrix] [m] (fmap g/simplify m))
#object[clojure.lang.MultiFn 0x41843e3e "
clojure.lang.MultiFn@41843e3e"
]
(defmethod g/invert [::matrix] [m] (invert m))
#object[clojure.lang.MultiFn 0x165df7c1 "
clojure.lang.MultiFn@165df7c1"
]
(defmethod g/make-rectangular [::matrix ::matrix] [a b]
(elementwise g/make-rectangular a b))
#object[clojure.lang.MultiFn 0x41c0c051 "
clojure.lang.MultiFn@41c0c051"
]
(defmethod g/make-polar [::matrix ::matrix] [a b]
(elementwise g/make-polar a b))
#object[clojure.lang.MultiFn 0x73349281 "
clojure.lang.MultiFn@73349281"
]
(defmethod g/real-part [::matrix] [m] (fmap g/real-part m))
#object[clojure.lang.MultiFn 0x4ed57a82 "
clojure.lang.MultiFn@4ed57a82"
]
(defmethod g/imag-part [::matrix] [m] (fmap g/imag-part m))
#object[clojure.lang.MultiFn 0x41ca178e "
clojure.lang.MultiFn@41ca178e"
]
(defmethod g/conjugate [::matrix] [m] (fmap g/conjugate m))
#object[clojure.lang.MultiFn 0x7fd6127a "
clojure.lang.MultiFn@7fd6127a"
]
(defmethod g/transpose [::matrix] [m] (transpose m))
#object[clojure.lang.MultiFn 0x2ecf7620 "
clojure.lang.MultiFn@2ecf7620"
]
(defmethod g/determinant [::square-matrix] [m] (determinant m))
#object[clojure.lang.MultiFn 0x213217f1 "
clojure.lang.MultiFn@213217f1"
]
(defmethod g/determinant [::s/structure] [s]
(two-tensor-> s (fn [m _] (determinant m))))
#object[clojure.lang.MultiFn 0x213217f1 "
clojure.lang.MultiFn@213217f1"
]
(defmethod g/trace [::square-matrix] [m] (trace m))
#object[clojure.lang.MultiFn 0x3ca928c1 "
clojure.lang.MultiFn@3ca928c1"
]
(defmethod g/trace [::s/structure] [s]
(two-tensor-> s (fn [m _] (trace m))))
#object[clojure.lang.MultiFn 0x3ca928c1 "
clojure.lang.MultiFn@3ca928c1"
]
(defmethod g/invert [::s/structure] [a]
(s:invert a))
#object[clojure.lang.MultiFn 0x165df7c1 "
clojure.lang.MultiFn@165df7c1"
]
(defmethod g/div [::s/structure ::s/structure] [rv s]
(s:divide-by-structure rv s))
#object[clojure.lang.MultiFn 0x233dab3a "
clojure.lang.MultiFn@233dab3a"
]
(defmethod g/solve-linear [::square-matrix ::s/up] [A b] (rsolve b A))
#object[clojure.lang.MultiFn 0x48e9c356 "
clojure.lang.MultiFn@48e9c356"
]
(defmethod g/solve-linear [::square-matrix ::s/down] [A b] (rsolve b A))
#object[clojure.lang.MultiFn 0x48e9c356 "
clojure.lang.MultiFn@48e9c356"
]
(defmethod g/solve-linear [::square-matrix ::column-matrix] [A b] (rsolve b A))
#object[clojure.lang.MultiFn 0x48e9c356 "
clojure.lang.MultiFn@48e9c356"
]
(defmethod g/solve-linear [::square-matrix ::row-matrix] [A b] (rsolve b A))
#object[clojure.lang.MultiFn 0x48e9c356 "
clojure.lang.MultiFn@48e9c356"
]
(defmethod g/solve-linear [::s/structure ::s/structure] [s product]
(s:solve-linear-left s product))
#object[clojure.lang.MultiFn 0x48e9c356 "
clojure.lang.MultiFn@48e9c356"
]
(defmethod g/solve-linear [::s/structure ::v/scalar] [s c]
(s/structure*scalar (s:invert s) c))
#object[clojure.lang.MultiFn 0x48e9c356 "
clojure.lang.MultiFn@48e9c356"
]
(defmethod g/solve-linear-right [::row-matrix ::square-matrix] [b A] (rsolve b A))
#object[clojure.lang.MultiFn 0x738d50cf "
clojure.lang.MultiFn@738d50cf"
]
(defmethod g/solve-linear-right [::down ::square-matrix] [b A] (rsolve b A))
#object[clojure.lang.MultiFn 0x738d50cf "
clojure.lang.MultiFn@738d50cf"
]
(defmethod g/solve-linear-right [::s/structure ::s/structure] [product s]
(s:solve-linear-right product s))
#object[clojure.lang.MultiFn 0x738d50cf "
clojure.lang.MultiFn@738d50cf"
]
(defmethod g/solve-linear-right [::v/scalar ::s/structure] [c s]
(s/scalar*structure c (s:invert s)))
#object[clojure.lang.MultiFn 0x738d50cf "
clojure.lang.MultiFn@738d50cf"
]
(defmethod g/dimension [::square-matrix] [m] (dimension m))
#object[clojure.lang.MultiFn 0x56adc8b2 "
clojure.lang.MultiFn@56adc8b2"
]
(defmethod g/dimension [::column-matrix] [m] (num-rows m))
#object[clojure.lang.MultiFn 0x56adc8b2 "
clojure.lang.MultiFn@56adc8b2"
]
(defmethod g/dimension [::row-matrix] [m] (num-cols m))
#object[clojure.lang.MultiFn 0x56adc8b2 "
clojure.lang.MultiFn@56adc8b2"
]

Column / Row Matrices only...

(defmethod g/dot-product [::row-matrix ::row-matrix] [a b]
(g/dot-product (row-matrix->down a)
(row-matrix->down b)))
#object[clojure.lang.MultiFn 0x4f7d2d76 "
clojure.lang.MultiFn@4f7d2d76"
]
(defmethod g/dot-product [::column-matrix ::column-matrix] [a b]
(g/dot-product (column-matrix->up a)
(column-matrix->up b)))
#object[clojure.lang.MultiFn 0x4f7d2d76 "
clojure.lang.MultiFn@4f7d2d76"
]
(defmethod g/inner-product [::row-matrix ::row-matrix] [a b]
(g/inner-product (row-matrix->vector a)
(row-matrix->vector b)))
#object[clojure.lang.MultiFn 0x31b14dc1 "
clojure.lang.MultiFn@31b14dc1"
]
(defmethod g/inner-product [::column-matrix ::column-matrix] [a b]
(g/inner-product (column-matrix->up a)
(column-matrix->up b)))
#object[clojure.lang.MultiFn 0x31b14dc1 "
clojure.lang.MultiFn@31b14dc1"
]
(defmethod g/cross-product [::row-matrix ::row-matrix] [a b]
(by-rows
(s/structure->vector
(g/cross-product (row-matrix->vector a)
(row-matrix->vector b)))))
#object[clojure.lang.MultiFn 0x7f10d912 "
clojure.lang.MultiFn@7f10d912"
]
(defmethod g/cross-product [::column-matrix ::column-matrix] [a b]
(up->column-matrix
(g/cross-product (column-matrix->up a)
(column-matrix->up b))))
#object[clojure.lang.MultiFn 0x7f10d912 "
clojure.lang.MultiFn@7f10d912"
]
(defmethod g/outer-product [::column-matrix ::row-matrix] [a b] (mul a b))
#object[clojure.lang.MultiFn 0x4ab81d57 "
clojure.lang.MultiFn@4ab81d57"
]
(defmethod g/partial-derivative [::matrix v/seqtype] [M selectors]
(fmap #(g/partial-derivative % selectors)
M))
#object[clojure.lang.MultiFn 0x3d9e83a6 "
clojure.lang.MultiFn@3d9e83a6"
]