neaderthal-linalg.txt

(ns mlstudy.linalgbook
  (:require [uncomplicate.neanderthal.linalg :as linalg]
            [tech.viz.vega :as vega])
  (:use     [uncomplicate.neanderthal core native vect-math]
            [clojure.repl]))

;;; See also
;; https://dragan.rocks/articles/21/Hello-World-Programming-with-Linear-Algebra
;; Based on Gareth Williams 9th Edition - Linear Algebra with Applications 2017

;; add up elements of a vector
(sum (dv [1 2 3 -1]))
;; => 5.0

;; Pg. 28 of LAA
;; Let u = -1,4,3,7 and v = -2,-3,1,0 be vectors inR4
;; Find u + v and 3u

;; u + v
;; sum vectors "x pLUS y"
(xpy (dv [-1 4 3 7]) (dv [-2 -3 1 0]))
;;    [  -3.00    1.00    4.00    7.00 ]

;; scalar multiplication "a times x"  - also "scal"
(ax 3 (dv [-1 4 3 7]))
;;    [  -3.00   12.00    9.00   21.00 ]

;; vector subtraction = u + -v => -1 * v + u = u -v
(axpy -1 (dv [1 2 3 4]) (dv [4 4 4 4]))
;; => #RealBlockVector[double, n:4, offset: 0, stride:1]
;;    [   3.00    2.00    1.00    0.00 ]

;; convenience x minus y
(defn my-xmy [x y] "x - y" (axpy -1 y x))

(my-xmy (dv [4 4 4 4]) (dv [1 2 3 4]))
;; => #RealBlockVector[double, n:4, offset: 0, stride:1]
;;    [   3.00    2.00    1.00    0.00 ]

;; pg 31. linear combination
;; u = [2 5 -3] v = [-4 1 9] w = [4 0 2]
;; what is 2u - 3v +w ?

(let [u (dv [2 5 -3])
      v (dv [-4 1 9])
      w (dv [4 0 2])]
  (xpy (ax 2 u) (ax -3 v) w))
;;    [  20.00    7.00  -31.00 ]


;; pg 32. Determine whether the vector 8, 0, 5 is a linear combination
;; of the vectors [1, 2, 3] [0, 1, 4] and [2, -1, 1]
;; answer: yes w/ 2 -1 3

;; coefficient matrix
@(def cma1 (dge 3 3 [1 2 3
                     0 1 4
                     2 -1 1]))
;;       →       1.00    0.00    2.00         
;;       →       2.00    1.00   -1.00         
;;       →       3.00    4.00    1.00         

;; right sides matrix?
@(def rsmb1 (dge 3 1 [8 0 5]))
;;       →       8.00         
;;       →       0.00         
;;       →       5.00         

;; sv a b: solve a system of linear equations
;; with coefficient matrix a and right sides matrix b
;; note sv already defined in native
(linalg/sv cma1 rsmb1)
;;       →       2.00         
;;       →      -1.00         
;;       →       3.00         

;; CORRECTAMUNDO!!!

;; Neater
(let [cma1 (dge 3 3 [1 2 3
                     0 1 4
                     2 -1 1])
      rsmb1 (dge 3 1 [8 0 5])]
  (linalg/sv cma1 rsmb1))
;;       →       2.00         
;;       →      -1.00         
;;       →       3.00         

;; Note - this looks similar to the above
;; From https://dragan.rocks/articles/17/Clojure-Linear-Algebra-Refresher-Linear-Transformations
(let [a (dge 3 3 [1 0 1 1 -1 1 3 1 4])
      b (dge 3 1 [11 -2 9])]
  (linalg/sv! a b)
  (def a-solution (col b 0))
  a-solution)
;;    [  17.00   -0.00   -2.00 ]


;; Pg. 45 Dot Product  1*1 + 2*2 + 3*3 = 1+4+9 = 14
(dot (dv 1 2 3) (dv 1 2 3))
;; => 14.0

;; commutative => order doesn't matter
;; associative => (4+5)+6 = 4+(5+6)
;; distributive => 6(5-2) = 6*5-6*2  i.e. 6*3 = 30-12



;; P 46. Norm (length) = sqrt (a^2 + b^2 .. n^2)

;; sqrt(4+9)=sqrt (13)
(nrm2 (dv 2 3)) 
;; => 3.605551275463989

(nrm2 (dv 3 4))
;; => 5.0

;; norm also = sqrt dot of u.u
(Math/sqrt (dot (dv 3 4) (dv 3 4)))
;; => 5.0

;; P. 47 Angle between vectors

;; cos(theta) = u.v / (norm u)*(norm v)
;; Determine the angle between the vectors u 1, 0, 0 and v 1, 0, 1
;; Answer: pi/4 = 0.7853981633974483  (45 degrees)

(let [u (dv 1 0 0)
      v (dv 1 0 1)]
  (Math/acos (/ (dot u v) (* (nrm2 u) (nrm2 v)))))
;; => 0.7853981633974484

;; so as a function

(defn angle
  "Angle between vectors u and v"
  [u v]
  (Math/acos (/ (dot u v) (* (nrm2 u) (nrm2 v)))))

(angle (dv 1 0 0) (dv 1 0 1))
;; => 0.7853981633974484

;; pg. 41 u and v are orthogonal if u.v = 0

;; Are (dv 2 -3 1) and (dv 1 2 4) orthogonal?
(zero? (dot (dv 2 -3 1) (dv 1 2 4)))
;; => true

;; Are (dv 1 2 -5) and (dv -1 3 1) orthogonal?
(zero? (dot (dv 1 2 -5) (dv -1 3 1)))
;; => true

;; orthonormal set : all pairs in set are orthogonal
;; standard basis for Rn is an orthonormal set
;; e.g. (1 0 0) (0 1 0) (0 0 1)

;; pg. 51 distance between points x & y
;; sqrt( (x1-y1)^2 + (x2-y2)^2 .. (xn-yn)^2)

;; example: what is the distance between (dv 1 -2 3 0) and (dv 4 0 -3 5)
;; answer sqrt(74) = 8.602325267042627

;; use x minus y
@(def vdf (my-xmy (dv 1 -2 3 0) (dv 4 0 -3 5)))
;;    [  -3.00   -2.00    6.00   -5.00 ]

(Math/sqrt (apply + (map #(* % %) vdf)))
;; => 8.602325267042627

;; in a func
(defn distance
  "Distance between vectors"
  [x y]
  (Math/sqrt (apply + (map #(* % %) (my-xmy x y)))))

(distance (dv 1 -2 3 0) (dv 4 0 -3 5))
;; => 8.602325267042627

;; show distance x,y = distance y,x
(distance (dv 4 0 -3 5) (dv 1 -2 3 0))
;; => 8.602325267042627

;; distance to self = 0
(distance (dv 4 0 -3 5) (dv 4 0 -3 5))
;; => 0.0


;; p. 58 curve fitting

;; example 1: polynomial degree 2 through points 1,6 2,3 3,2
;; polynomial is y = a0 + a1x + a2x^2

;; coefficients
(dge 3 3 [1 1 1 1 2 4 1 3 9] {:layout :row})
;;       →       1.00    1.00    1.00         
;;       →       1.00    2.00    4.00         
;;       →       1.00    3.00    9.00         

(let [cma1 (dge 3 3 [1 1 1
                     1 2 4
                     1 3 9] {:layout :row})
      rsmb1 (dge 3 1 [6 3 2])]
  (linalg/sv cma1 rsmb1))
;;       →      11.00         
;;       →      -6.00         
;;       →       1.00         

;; result checks out
;; dge defaul layout is column. calling cols gets the right order
(map (partial dot (dv 11 -6 1)) (cols (dge 3 3 [1 1 1 1 2 4 1 3 9])))
;; => (6.0 3.0 2.0)

;; so parabola through the points is y = 11 -6x + x^2


;; p. 59 Electrical Network Analysis

(let [cma1 (dge 3 3 [1 1 -1
                     4 0 1
                     0 4 1] {:layout :row})
      rsmb1 (dge 3 1 [0 8 16])]
  (linalg/sv cma1 rsmb1))
;;       →       1.00         
;;       →       3.00         
;;       →       4.00         

;; currents are I1 = 1, I2 = 3, I3 = 4 amps

;; example 3
(let [cma1 (dge 3 3 [1 1 -1
                     1 0 2
                     1 -2 0] {:layout :row})
      rsmb1 (dge 3 1 [0 12 -4])]
  (linalg/sv cma1 rsmb1))
;;       →       2.00         
;;       →       3.00         
;;       →       5.00         

;; currents are I1 = 2, I2 = 3, I3 = 5 amps


;;;; CHAPTER TWO

;; p. 70 Matrices and Linear Transformations

(= (dge 3 3 [1 2 3 4 5 6 7 8 9]) (dge 3 3 [1 2 3 4 5 6 7 8 9]))
;; => true

;; add matrices
(xpy (dge 2 3 [1 0 4 -2 7 3]) (dge 2 3 [2 -3 5 1 -6 8]))
;;       →       3.00    9.00    1.00         
;;       →      -3.00   -1.00   11.00         

;; scalar multiplication

;; note vector is col-1 then col-2 etc.
(dge 2 3 [1 7 -2 -3 4 0])
;;       →       1.00   -2.00    4.00         
;;       →       7.00   -3.00    0.00         

(ax 3 (dge 2 3 [1 7 -2 -3 4 0]))
;;       →       3.00   -6.00   12.00         
;;       →      21.00   -9.00    0.00         

;; negation
(ax -1 (dge 2 3 [1 -3 0 6 -7 2]))
;;       →      -1.00    0.00    7.00         
;;       →       3.00   -6.00   -2.00         

;; subtraction
(my-xmy (dge 2 3 [5 3 0 6 -2 -5]) (dge 2 3 [2 0 8 4 -1 6]))
;;       →       3.00   -8.00   -1.00         
;;       →       3.00    2.00  -11.00         

;; p. 72 Matrix Multiplication

;; example 1: A:2x2 B:2x3 results in 2x3
(mm (dge 2 2 [1 2 3 0]) (dge 2 3 [5 3 0 -2 1 6]))
;;       →      14.00   -6.00   19.00         
;;       →      10.00    0.00    2.00         

;; size of product matrix
;; AB exists of size m x n for matrices A: m x r and B: r x n

;; my example : A:2x3 B:3x4 results in 2x4
(mm (dge 2 3 [1 2 3 4 5 0]) (dge 3 4 [5 3 0 -2 1 6 1 2 3 4 5 6]))
;;       →      14.00   31.00   22.00   49.00         
;;       →      22.00    0.00   10.00   28.00         

;; p. 74 special matrices

;; zero matrix - just don't provide values to dge
(dge 3 5)
;;       →       0.00    0.00    0.00    0.00    0.00         
;;       →       0.00    0.00    0.00    0.00    0.00         
;;       →       0.00    0.00    0.00    0.00    0.00         

;; submatrix
(def subex (dge 3 4 [5 3 0 -2 1 6 1 2 3 4 5 6]))
;;       →       5.00   -2.00    1.00    4.00         
;;       →       3.00    1.00    2.00    5.00         
;;       →       0.00    6.00    3.00    6.00         

;; submatrix returns the submatrix of matrix `a` starting from row `i`, column `j`,
;; that has `k` rows and `l` columns.

;; submatrix starting at 0,1 of size 2x2 
(submatrix subex 0 1 2 2)
;;       →      -2.00    1.00         
;;       →       1.00    2.00         

;; matrix algebraic properties
;; A + B = B + A : relexive
;; A + (B + C) = (A + B) + C : associative under addition
;; additive identity 0
;; distributive under addition

;; A (BC) = (AB)C : associative under multiplication
;; distributive under multiplication
;; mutiplicative identity 1
;; scalar multiplcation is also associaaative and distributive

;; NOTE!!
;; AB != BA in general!
;; AB = AC does not imply that B = C.
;; PQ = 0 does not imply that P = O or Q = O.

;; powers of matrices

;; A^k = AAA..A k times. A^rA^s = A^(r+s)
;; (A^r)^s = A^(r*s)
;; A^0 = Isubn - identity (by definition)

;; p. 92 Transpose of a matrix

;; example 1 : tranpose [2 -8 7 0]
@(def amat1 (dge 2 2 [2 -8 7 0]))
;;       →       2.00    7.00         
;;       →      -8.00    0.00         

(trans amat1)
;;       →       2.00   -8.00         
;;       →       7.00    0.00         

;; another
@(def amat2 (dge 2 3 [1 4 2 5 -7 6]))
;;       →       1.00    2.00   -7.00         
;;       →       4.00    5.00    6.00         

(trans amat2)
;;       →       1.00    4.00         
;;       →       2.00    5.00         
;;       →      -7.00    6.00         

;; definition: symmetric matrix is equal to its transpose
@(def symmat (dge 3 3 [0 1 -4 1 7 8 -4 8 3]))
;;       →       0.00    1.00   -4.00         
;;       →       1.00    7.00    8.00         
;;       →      -4.00    8.00    3.00         

(= symmat (trans symmat))
;; => true

;; Let A and B be symmetric matrices of the same size.
;; The product AB is symmetric if and only if AB = BA.

;; p. 95 trace of a (square) matrix is the sum of its diagonal elements

@(def sqm (dge 3 3 [4 2 7 1 -5 3 -2 6 0]))
;;       →       4.00    1.00   -2.00         
;;       →       2.00   -5.00    6.00         
;;       →       7.00    3.00    0.00         

;; no "trace" func
(dia sqm)
;;    [   4.00   -5.00    0.00 ]

(sum (dia sqm))
;; => -1.0

;; so convenience func trace
(defn trace
  "Return the trace (sum of diagonal elements)  of m"
  [m]
  (sum (dia m)))

(trace sqm)
;; => -1.0

;; p. 102 Section 2.4 : Inverse of a Matrix

;; A is invertible if AB = BA = I for some matrix B
;; prove [1 3 2 4] has inverse [-2 3/2 1 -1/2]
@(def invm (dge 2 2 [1 3 2 4]))
;;       →       1.00    2.00         
;;       →       3.00    4.00         

@(def invmi (dge 2 2 [-2 3/2 1 -1/2]))
;; => #RealGEMatrix[double, mxn:2x2, layout:column, offset:0]
;;       ▥       ↓       ↓       ┓    
;;       →      -2.00    1.00         
;;       →       1.50   -0.50         
;;       ┗                       ┛    

;; trf Triangularizes a non-TR matrix
;; tri Computes the inverse of a triangularized matrix
@(def invmi (linalg/tri (linalg/trf invm)))
;;       →      -2.00    1.00         
;;       →       1.50   -0.50         

;; check that AAinv = I (it does)
(mm invm invmi)
;;       →       1.00    0.00         
;;       →       0.00    1.00         

;; Theorem 2.7 inverses are unique
;; p. 106
;; Theorem 2.8
;; Let AX = Y be a system of n linear equations in n variables.
;; If Ainv exists, the solution is unique and is given by X = AinvY.

;; cryptography

;; Encryption Matrix
@(def encryption-matrix (dge 3 3 [-3 0 4
                                  -3 1 3
                                  -4 1 4]))
;;       →      -3.00   -3.00   -4.00         
;;       →       0.00    1.00    1.00         
;;       →       4.00    3.00    4.00         

;; Message (as a matrix)
;;                  b u  y  -  i b m  -  s  t  o  c k  -  -
@(def msg (dge 3 5 (map int "buy ibm stock  ")))
;;       →       2.00   27.00   13.00   20.00   11.00         
;;       →      21.00    9.00   27.00   15.00   27.00         
;;       →      25.00    2.00   19.00    3.00   27.00         

@(def ciphtxt (mm encryption-matrix msg))
;;       →    -169.00 -116.00 -196.00 -117.00 -222.00         
;;       →      46.00   11.00   46.00   18.00   54.00         
;;       →     171.00  143.00  209.00  137.00  233.00         

;; to decrypt we need the inverse of the encryption matrix
@(def decryption-matrix (linalg/tri (linalg/trf encryption-matrix)))
;;       →       1.00    0.00    1.00         
;;       →       4.00    4.00    3.00         
;;       →      -4.00   -3.00   -3.00         

@(def clrtxt (mm decryption-matrix ciphtxt))
;;       →       2.00   27.00   13.00   20.00   11.00         
;;       →      21.00    9.00   27.00   15.00   27.00         
;;       →      25.00    2.00   19.00    3.00   27.00         

;; was decryption successful?
(= msg clrtxt)
;; => true

;; show it
(apply str (map char (apply concat clrtxt)))
;; => "buy ibm stock  "


;; p. 116 Transformations

;; scaling / dilation and contraction

;; scale by 2 using matrix (dilation)
(mv (dge 2 2 [2 0 0 2]) (dv 3 4))
;;    [   6.00    8.00 ]

;; scale by -1/2 (contraction)
(mv (dge 2 2 [-1/2 0 0 -1/2]) (dv 3 4))
;;    [  -1.50   -2.00 ]

;; these are equivalent to scalar multiply, e.g.
(ax 2 (dv [3 4]))
;;    [   6.00    8.00 ]

;; p. 117 Reflection - maps every point to its 'mirror'
(mv (dge 2 2 [1 0 0 -1]) (dv 3 4))
;;    [   3.00   -4.00 ]


;; p. 118 rotation about the origin e.g. 180 degrees
;; convenience
(defn _sin [th] (Math/sin th))
(defn _-sin [th] (- (Math/sin th)))
(defn _cos [th] (Math/cos th))
(defn _-cos [th] (- (Math/cos th)))

(def theta (/ Math/PI 2)) 

(mv (dge 2 2 [(_cos theta) (_sin theta)
              (_-sin theta) (_cos theta)]) (dv 3 2))
;;    [  -2.00    3.00 ]

;; let's make a rotation function
(defn rotv [v theta]
  (mv (dge 2 2 [(_cos theta) (_sin theta)
              (_-sin theta) (_cos theta)]) (dv v)))

;; 90 degrees
(rotv [3 2] theta)
;;    [  -2.00    3.00 ]

;; 45 degrees
(rotv [3 2] (/ Math/PI 4))
;;    [   0.71    3.54 ]

;; also 45 degrees
(rotv [1 0] (/ Math/PI 4))
;; => #RealBlockVector[double, n:2, offset: 0, stride:1]
;;    [   0.71    0.71 ]

;; also 45 degrees
(rotv [(Math/sqrt 2) 0] (/ Math/PI 4))
;; => #RealBlockVector[double, n:2, offset: 0, stride:1]
;;    [   1.00    1.00 ]

;; MATRIX TRANSFORMATIONS

;;  Matrix transformations map line segments into line segments (or
;;  points). If the matrix is invertible, the transformation also maps
;;  parallel lines into parallel lines.

;; p. 120 Composition of Transformations

;; The composite transformation T = Tn...T1 is defined by the matrix
;; product An...A1 (assuming this product exists).

;; composite of: reflection, rotation, dilation
@(def comprotmat 
  (let [reflmat (dge 2 2 [3 0 0 3])
        rotmat  (dge 2 2 [(_cos (/ Math/PI 2))
                          (_sin (/ Math/PI 2))
                          (_-sin (/ Math/PI 2))
                          (_cos (/ Math/PI 2))])
        dilmat (dge 2 2 [1 0 0 -1])]
    (mm reflmat rotmat dilmat)))
;;       →       0.00    3.00         
;;       →       3.00   -0.00         

(mv comprotmat (dv 1 2))
;;    [   6.00    3.00 ]


;; p. 121 Theorem 2.9 Orthogonal Transformations.

;; Orthogonal transformations preserve norms, angles, and distances.
;; Orthogonal transformations preserve the shapes of rigid bodies and
;; are often referred to as rigid motions.

;; p. 123 Translation

;; translation is just vector addition
(xpy (dv [1 2]) (dv [4 2]))
;;    [   5.00    4.00 ]

;; p. 124 Affine Transformation
;; T(u) = A(u) + v

(defn affinexform [am av v]
  (xpy
   (mv (dge 2 2 am) (dv v))
   (dv av)))

(def myaffine (partial affinexform [2 1 1 1] [1 2]))

(myaffine [1 0])
;; => #RealBlockVector[double, n:2, offset: 0, stride:1]
;;    [   3.00    3.00 ]

(def unitsq [[1 0] [1 1] [0 1] [0 0]])
(map myaffine unitsq)
;;    [   3.00    3.00 ]
;;    [   4.00    4.00 ]
;;    [   2.00    3.00 ]
;;    [   1.00    2.00 ]


;; p. 126 Linear Transformations, Graphics, and Fractals

;; Linear Transformations preserve addition and scalar multiplication
;; Every matrix transformation is linear (dilation, contraction,
;; reflection, rotation)

;; p. 132 Example 5 Rotate about a point
;; use 'homogenous coordinates' which are 3x3 tranformation matrices
;; translate rotation point to the origin, rotate, translate back.

(defn rotapointmx [h k theta]
  (let [txorig (dge 3 3 [1 0 0 0 1 0 (- h) (- k) 1])
        rot (dge 3 3 [(_cos theta) (_sin theta) 0
                      (_-sin theta) (_cos theta) 0
                      0 0 1])
        txback (dge 3 3 [1 0 0 0 1 0 h k 1])]
    (mm txback rot txorig)))

@(def myrotapt (rotapointmx 5 4 (/ Math/PI 2)))
;; => #RealGEMatrix[double, mxn:3x3, layout:column, offset:0]
;;       ▥       ↓       ↓       ↓       ┓    
;;       →       0.00   -1.00    9.00         
;;       →       1.00    0.00   -1.00         
;;       →       0.00    0.00    1.00         
;;       ┗                               ┛    


;; p. 134 Fractals

;; see ~/gitrepos/clojure/fracleaf
;; the pct rand function call stuff has been extracted to randfunpct.clj


;; p. 143: A stochastic matrix is a square matrix whose elements are
;;  probabilities and whose columns add up to 1.

;; THEOREM 2.10: If A and B are stochastic matrices of the same size,
;; then AB is a stochastic matrix.
 
;; p 145 Eample 1 - land use trends

;; This type of model, where the probability of going from one state
;; to another depends only on the current state rather than on a more
;; complete historical description, is called a Markov Chain.*

;; matrix of transition probabilities

;; city - suburb transition matrix
(def csub-txn (dge 2 2 [0.97 0.03 0.01 0.99]))

;; starting population [city, country]
(def stpop (dv 83 178))

;; after one year, [city, country]
(mv csub-txn stpop)
;; => #RealBlockVector[double, n:2, offset: 0, stride:1]
;;    [  82.29  178.71 ]

(defn nxtyrpop [yr] (mv csub yr))

;; 5 years of migrations between cities and suburbs
(take 5 (iterate nxtyrpop stpop))
;;    [  83.00  178.00 ]
;;    [  82.29  178.71 ]
;;    [  81.61  179.39 ]
;;    [  80.95  180.05 ]
;;    [  80.33  180.67 ]

;; p 147 Guinea Pigs

;; guinea pig generation transition matrix
@(def gp-generation (dge 3 3 [1 0 0 1/2 1/2 0 0 1 0]))
@(def gp-initial (dv 1/3 1/3 1/3))

(defn nxtgengp [gen] (mv gp-generation gen))

;; 5 generations of guinea pigs
(take 5 (iterate nxtgengp gp-initial))
;;    [   0.33    0.33    0.33 ]
;;    [   0.50    0.50    0.00 ]
;;    [   0.75    0.25    0.00 ]
;;    [   0.88    0.13    0.00 ]
;;    [   0.94    0.06    0.00 ]

;; p. 149 A Communication Model

;; A digraph is a finite collection of vertices P1, P2, c, Pn,
;; together with directed arcs joining cer- tain pairs of vertices. A
;; path between vertices is a sequence of arcs that allows one to
;; proceed in a continuous manner from one vertex to another. The
;; length of a path is its number of arcs. A path of length n is
;; called an n-path.

;; A digraph can be described by a matrix A, called its adjacency
;; matrix. This matrix consists of zeros and ones and is defined as
;; follows.

;; a(ij) = 1 if arc from p(i) to p(j)
;;       = 0 if no arc from p(i) to p(j)
;;       = 0 if i = j

;; THEOREM 2.11 If A is the adjacency matrix of a digraph, let aij^(m) be the element in row i and column j of A^(m). The number of m-paths from Pi to Pj = a(ij)^(m).

;; p 151 powers of adjacency matrix

(def amatx (dge 5 5 [0 1 0 0 0
                     1 0 1 0 1
                     0 0 0 1 0
                     0 0 0 0 1
                     0 1 1 1 0]))

;; how many 2 paths from i to j = a2(i,j)
;; e.g. zero 2 paths from 2 to 3, 2 2 paths from 4 to 2
@(def a2 (mm amatx amatx))
;;       →       1.00    0.00    0.00    0.00    1.00         
;;       →       0.00    2.00    0.00    1.00    0.00         
;;       →       1.00    1.00    0.00    1.00    1.00         
;;       →       0.00    2.00    0.00    1.00    1.00         
;;       →       1.00    0.00    1.00    0.00    2.00         

;; how many 3 paths from i to j = a3(i,j)
@(def a3 (mm amatx a2))
;;       →       0.00    2.00    0.00    1.00    0.00         
;;       →       2.00    0.00    1.00    0.00    3.00         
;;       →       1.00    2.00    1.00    1.00    2.00         
;;       →       2.00    1.00    1.00    1.00    3.00         
;;       →       0.00    4.00    0.00    2.00    1.00         

;; etc.

;; This model that we have discussed gives the lengths of the shortest
;; paths of a digraph; it does not give the intermediate stations that
;; make up that path. Mathematicians have not, as yet, been able to
;; derive this information from the adjacency matrix. An algorithm for
;; finding the shortest paths for a specific digraph, using a search
;; procedure, has been developed by a Dutch computer scientist named
;; Edsger Dijkstra.

;; p 152 - Distance in a Digraph

;; d(ij) = # arcs in shortest path from i to j
;;       = 0 if i = j
;;       = x(undefined) if there is no path from i to j


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;; CHAPTER 3 - Determinants and Eigenvectors ;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Every square matrix has a 'determinant'

;; Eigenvalues and eigenvectors are special scalars and vectors
;; associated with matri- ces. They can be expressed in terms of
;; determinants.

;; 2x2: determinant of a =  a11*a22 - a12*a21

@(def trri (linalg/trf (dge 2 2 [2 -3 4 1])))
;;       →      -3.00    1.00         
;;       →      -0.67    4.67         

(linalg/det trri)
;; => -14.0

(dge 2 2 [2 -3 4 1])
;;       →       2.00    4.00         
;;       →      -3.00    1.00         

(->> (dge 2 2 [2 -3 4 1])
     linalg/trf
     linalg/det)
;; => -14.0

;; That's wrong - should be 14:
;; 2*1 - -3*4 = 2 - -12 = 14

;;; NOTED BUG REPORTED and FIXED!

;; determinants "give us an algebraic formula for the inverse of an arbitrary matrix"

;; the determinant of a 2 3 2 matrix is given by the difference of the products of the two diagonals of the matrix.

;; p. 165 The determinant of a square matrix is the sum of the
;; products of the elements of the first row and their cofactors

;; try w/ 3x3
@(def trri2 (linalg/trf (dge 3 3 [1 3 4 2 0 2 -1 1 1])))
;;       →       4.00    2.00    1.00         
;;       →       0.75   -1.50    0.25         
;;       →       0.25   -1.00   -1.00         

(linalg/det trri2)
;; => -6.0
;; correct

@(def trri3 (linalg/trf (dge 4 4 [2  0  7  0
                                  1 -1 -2  1
                                  0  0  3  0
                                  4  2  5 -3])))
;;       →       7.00   -2.00    3.00    5.00         
;;       →       0.29    1.57   -0.86    2.57         
;;       →       0.00   -0.64   -0.55    3.64         
;;       →       0.00    0.64   -1.00   -1.00         

(linalg/det trri3)
;; => 6.000000000000003
;; 6 is correct (pretty close)

;; A square matrix A is said to be singular if det(A) = 0. A is nonsingular otherwise.
 
@(def trri4 (linalg/trf (dge 3 3 [3 -1 2
                                  4 -6 9
                                  -2 3 -3])))
;;       ▥       ↓       ↓       ↓       ┓    
;;       →       3.00    4.00   -2.00         
;;       →       0.67    6.33   -1.67         
;;       →      -0.33   -0.74    1.11         
(linalg/det trri4)
;; => -21.000000000000004

@(def trri5 (linalg/trf (dge 3 3 [1 0 -2
                                  4 2 -4
                                  3 5 10])))
;;       →      -2.00   -4.00   10.00         
;;       →      -0.00    2.00    5.00         
;;       →      -0.50    1.00    3.00         

(linalg/det trri5)
;; => 12.0


;; p 171 properties of determinants


;; Theorem 3.2

;; (a) if you multiply a rows(column) by c: |B| = c|A|
;; e.g. multiply row 1 x 2
(->> (dge 3 3 [2 0 -2
               8 2 -4
               6 5 10])
     linalg/trf
     linalg/det)
;; => 24.0 (confirmed)

;; (b) if you interchange rows: |B| = -|A|
(->> (dge 3 3 [0 1 -2
               2 4 -4
               5 3 10])
     linalg/trf
     linalg/det)
;; => -12.0 - confirmed (see trri5)

;; (c) add a multiple of one row(column) to another |B| = |A| (no change)
(->> (dge 3 3 [1 2 -2
               4 10 -4
               3 11 10])
     linalg/trf
     linalg/det)
;; => 11.999999999999996 (close!)


;; THEOREM 3.3
;; Let A be a square matrix. A is singular if
;; (a) all the elements of a row(column) are zero.
;; (b) two rows(columns) are equal.
;; (c) two rows(columns) are proportional.

;; THEOREM 3.4
;; Let A and B be nxn matrices and c be a nonzero scalar.
;; (a) Determinant of a scalar multiple: |cA| = C^n|A|
;; (b) Determinant of a product: |AB| = |A||B|
;; (c) Determinant of a transpose: |At| = |A|
;; (d) Determinant of an inverse: A^-1 = 1/|A|


;; p. 175 Upper Triangular

;; matrix is called an upper triangular matrix when all nonzero
;; elements lie on or above the main diagonal. Likewise with lower triangular

;; The determinant of a triangular matrix is the product of its diagonal elements.

;; p188 3.4 Eigenvalues and Eigenvectors

;; Let A be an nXn matrix. A scalar L is called an eigenvalue of A if there exists a nonzero vector x in Rn such that
;;                                 Ax = Lx
;; The vector x is called an eigenvector corresponding to L.

;; example
(def A (dge 2 2 [-4 3
                 -6 5]))

;; solution 1: L = 2
(def L1 2)
(def x1 (dv [-1 1]))  
(mv A x1)
;; => #RealBlockVector[double, n:2, offset: 0, stride:1]
;;    [  -2.00    2.00 ]

(ax L1 x1)
;; => #RealBlockVector[double, n:2, offset: 0, stride:1]
;;    [  -2.00    2.00 ]

;; confirmed

;; solution 2: L = -1
(def L2 -1)
(def x2 (dv [-2 1]))  

;; x is eigenvector if Ax = Lx
(mv A x2)
;; => #RealBlockVector[double, n:2, offset: 0, stride:1]
;;    [   2.00   -1.00 ]

(ax L2 x2)
;; => #RealBlockVector[double, n:2, offset: 0, stride:1]
;;    [   2.00   -1.00 ]

;; confirmed