Little matrix DSL

matrix.red

Red [
	Author: "Toomas Vooglaid"
	Date: 7-9-2017
	Last-update: 4-10-2017
]
mx: context [
	ctx: self
	mtx: object [
		rows: cols: data: none
		get-col: func [col][extract at data col cols]
		get-row: func [row][copy/part at data row - 1 * cols + 1 cols]
		remove-row: func [row][remove/part at data get-idx row 1 cols rows: rows - 1 show]
		remove-col: func [col][
			loop rows [data: remove skip data cols - 1]
			data: head data cols: cols - 1 show
		]
		get-idx: func [row col][index? at data row - 1 * cols + col]
		get-at: func [row col][pick data row - 1 * cols + col]
		get-row-idx: func [idx][idx - 1 / cols + 1]
		get-col-idx: func [idx][idx - 1 % cols + 1]
		get-diagonal: func [i dir /local out][
			data: skip data i - 1 
			set [comp inc] switch dir [r [0 :+] l [1 :-]]
			out: collect [
				while [not tail? data][
					keep data/1 
					data: case [
						all [dir = 'r 0 = ((index? data) % cols)] [next data]
						all [dir = 'l 1 = ((index? data) % cols)] [skip data 2 * cols - 1]
						true [skip data cols + either dir = 'r [1][-1]]
					]
				]
			]
			data: head data
			out
		]
		get-order: does [to-pair reduce [rows cols]]
		to-float: does [forall data [data/1: system/words/to-float data/1]]
		swap-dim: has [c][c: cols cols: rows rows: c]
		square?: does [rows = cols]
		symmetric?: has [d][transpose d: copy data transpose equal? data d]
		diagonal?: function [][
			either square? [
				repeat i cols [
					repeat j rows [
						if (i <> j) and (0 <> get-at i j) [return false]
				]] 
				true
			][false]
		]
		zero?: does [0 = ctx/summa data]
		singular?: degenerate?: does [determinant = 0]
		invertible?: nonsingular?: nondegenerate?: does [not singular?]
		sub-exclude: func [rs cs /local m2][ ; TBD
			m2: copy self
			switch type?/word rs [
				block! [
					sort/reverse rs
					forall rs [
						remove
					]
				]
				integer! []
				none! []
			]
			switch type?/word cs [
				block! []
				integer! []
				none! []
			]
		]
		transpose: does [ctx/transpose self]
		rotate: func [n][ctx/rotate n self]
		show: does [new-line/skip copy data true cols]
		pretty: function [][;/bar /local d i col-lengths][
			col-lengths: copy []
			repeat i cols [
				c: copy get-col i
				c: sort/compare c func [a b][(length? form a) > (length? form b)]
				append col-lengths length? form first c
			]
			cols2: copy []
			templ: copy []
			letters: "abcdefghijklmnopqrstuvwyz" 
			repeat n cols [
				append cols2 to-word pick letters n
				append templ compose [
					pad/left (to-word pick letters n) (pick col-lengths n) 
					;(either bar and (n < cols) ["│"][""])⎡⎢⎤⎣⎥⎦  ⎾⏋⎿⏌  ⌈⌉⌊⌋
				]
			]
			step: (summa col-lengths) + cols - 1 
			print [#"┌" pad #" " step #"┐"]
			foreach (cols2) data [
				print [#"│" reduce compose templ #"│"]
			]
			print [#"└" pad #" " step #"┘" #"^/"]
		]
		swap-rows: func [r1 r2][ctx/swap-rows r1 r2 self]
		determinant: does [ctx/determinant self]
		trace: does [ctx/trace self]
		identity: func [/side d][either side [ctx/identity/side self d][ctx/identity self]]
		split-col: func [col][ctx/split-col col self]
		neighbours: func [row col /local out][
			out: copy []
			foreach [r c] case [
				all [row = 1 col = 1][[1 2 2 1 2 2]]
				all [row = 1 col = cols] [reduce [
						1 cols - 1 2 cols - 1 2 cols
					]] 
				all [row = rows col = 1] [reduce [
						rows - 1 1 rows - 1 2 rows 2
					]]
				all [row = rows col = cols] [reduce [
						rows - 1 cols - 1 rows - 1 cols rows cols - 1
					]] 
				row = 1  [reduce [
						1 col - 1 1 col + 1 2 col - 1 2 col 2 col + 1
					]]
				row = rows [reduce [
						row - 1 col - 1 row - 1 col row - 1 col + 1 row col - 1 row col + 1
					]]
				col = 1 [reduce [
						row - 1 col row - 1 col + 1 row col + 1 row + 1 col row + 1 col + 1
					]]
				col = cols [reduce [
						row - 1 col - 1 row - 1 col row col - 1 row + 1 col - 1 row + 1 col
					]]
				true [reduce [
						row - 1 col - 1 row - 1 col row - 1 col + 1 row col - 1 row col + 1 row + 1 col - 1 row + 1 col row + 1 col + 1
					]]
			][append out get-at r c]
			out
		]
	]
	vector-op: func [op a b /local i][
		case [
			all [number? a number? b]	[return either op? :op [a op b][op a b]]
			all [number? a any-block? b][forall b [b/1: either op? :op [a op b/1][op a b/1]] return b]
			all [any-block? a number? b][forall a [a/1: either op? :op [a/1 op b][op a/1 b]] return a]
			all [any-block? a any-block? b][
				either (length? a) = (length? b) [
					forall a [i: index? a a/1: either op? :op [a/1 op pick b i][op a/1 pick b i]]
					return a
				][
					cause-error 'user 'message ["Vectors must be of the same length!"]
				]
			]
		]
	]
	product: func [blk /local out][out: 1 forall blk [out: out * blk/1] out]
	summa: func [blk /local out][out: 0 forall blk [out: out + blk/1] out]
	determinant: func [m /local i r l][
		either m/square? [
			switch/default m/cols [
				0	[1]
				1	[m/data/1]
				2 	[math [m/data/1 * m/data/4 - m/data/2 * m/data/3]]
			][
				r: make block! m/cols l: make block! m/cols
				repeat i m/cols [
					insert r product m/get-diagonal i 'r
					insert l product m/get-diagonal i 'l
				]
				(summa r) - (summa l)
			]
		][
			cause-error 'user 'message ["Matrix must be square to find determinant!"]
		]
	]
	trace: func [m][
		either m/square? [
			summa m/get-diagonal 1 'r
		][
			cause-error 'user 'message ["Trace is defined for square matrices only!"]
		]
	]
	add: func [op m1 m2][
		either all [m1/cols = m2/cols m1/rows = m2/rows][;length? m1/data length? m2/data [ 
			repeat i length? m1/data [m1/data/:i: m1/data/:i op m2/data/:i]
		][
			cause-error 'user 'message ["Matrices of unequal dimensions!"]
		]
		m1
	]
	multi: func [m1 m2 /local m3 val i j k l][
		either equal? l: m1/cols m2/rows [
			m3: make mtx [rows: m1/rows cols: m2/cols data: make block! (m1/rows * m2/cols)]
			repeat i m1/rows [
				repeat j m2/cols [
					val: 0
					repeat k l [val: (m1/get-at i k) * (m2/get-at k j) + val]
					append m3/data val
				]
			]
		][
			cause-error 'user 'message ["Dimensions don't match in multiplication!"]
		]
		m3
	]
	kronecker: func [m1 m2 /local m3 i j k l][
		m3: make mtx [rows: m1/rows * m2/rows cols: m1/cols * m2/cols data: make block! rows * cols]
		repeat i m1/rows [
			repeat j m2/rows [
				repeat k m1/cols [
					repeat l m2/cols [
						append m3/data (m1/get-at i k) * (m2/get-at j l)
		]]]]
		m3
	]
	transpose: func [m /local d i j r c][
		d: copy []
		repeat i c: m/cols [repeat j r: m/rows [append d m/get-at j i]]
		m/cols: r m/rows: c	m/data: d
		m
	]
	rotate: func [n [integer!] m /local data i][
		data: copy []
		switch n [
			1 or -3 [repeat i m/cols [append data copy reverse m/get-col i] m/swap-dim]
			2 or -2 [repeat i m/rows [append data reverse copy m/get-row m/rows + 1 - i]]
			3 or -1 [repeat i m/cols [append data copy m/get-col m/cols + 1 - i] m/swap-dim]
		]
		m/data: data 
		m
	]
	swap-rows: function [r1 r2 m][
		tmp: m/get-row r1 
		change/part at m/data r1 - 1 * m/cols + 1 m/get-row r2 m/cols 
		change/part at m/data r2 - 1 * m/cols + 1 tmp m/cols 
		m
	]
	identity: func [m /side d /local i][
		d: either side [case [d = 'l ['rows] d = 'r ['cols]]]['rows]
		either (side or m/square?) [
			data: make block! power m/:d 2
			repeat i m/:d [repeat j m/:d [append data either i = j [1][0]]]
			make mtx compose [cols: (m/:d) rows: (m/:d) data: (reduce [data])] 
		][
			cause-error 'user 'message ["You need to determine /side ['l | 'r] for non-square matrix!"]
		]
	]
	augment: func [m1 m2 /local i j][
		either m1/rows = m2/rows [
			repeat i m1/rows [
				k: m1/rows - i + 1
				j: m1/get-idx k m1/cols + 1
				insert at m1/data j m2/get-row k
			]
			m1/cols: m1/cols + m2/cols
		][
			cause-error 'user 'message ["Augmented matrix must have same number of rows as the other!"]
		]
		m1
	]
	rref: func [n /local m i j c val][
		m: copy/deep n m/to-float
		repeat i m/rows [
			; make the pivot
			if 0 = m/get-at i i [
				c: at m/get-col i i + 1
				until [
					c: next c 
					if tail? c [
						cause-error 'user 'message ["Impossible to get reduced row eschelon form!"]
					] 
					0 < first c
				]
				swap-rows i index? c 
			]
			; reduce it to 1
			if 1 <> val: m/get-at i i [
				change/part at m/data m/get-idx i 1 vector-op :/ m/get-row i val m/cols
			]
			; reduce other rows at this column to 0 
			repeat j m/rows [
				if all [j <> i 0 <> c: m/get-at j i][
					change/part at m/data m/get-idx j 1 vector-op :- m/get-row j vector-op :* c m/get-row i m/cols
				]
			]
		]
		m
	]
	split-col: func [col m /local data i j cls][
		data: copy []
		cls: m/cols - col + 1
		repeat i m/rows [
			j: m/rows - i + 1
			insert data take/part at m/data m/get-idx j col cls 
		] 
		m/cols: col - 1
		reduce [m make mtx compose/deep [rows: (m/rows) cols: (cls) data: [(data)]]]
	]
	invert: func [m /local n][
		augment m identity m
		n: rref m
		m: first split-col m/rows + 1 m
		second split-col n/rows + 1 n
	]
	game-of-life: func [m /local out c i][
		out: copy []
		c: m/data
		forall c [
			s: summa m/neighbours m/get-row-idx i: index? c m/get-col-idx i
			append out case [
				all [c/1 = 1 s < 2][0]
				all [c/1 = 1 find [2 3] s][1]
				all [c/1 = 1 s > 3][0]
				all [c/1 = 0 s = 3][1]
				true [0]
			]
		]
		m/data: out
		m
	]
	
	ops-rule: ['+ | '- | '* | '/ | '% | '** | '>> | '<< | '>>> | 'and | 'or | 'xor | 'div | 'x | 'augment]

	set 'matrix func [spec /local dim dims rule result m w m1 m2 op op' var vars ops unary unaries d matrices][
		vars: copy [] ops: copy [] matrices: copy [] unaries: copy []
		matrix-rule: [(m: none)[
			set dim pair! [set mdata block! | set w word! if (block? get/any w)(mdata: get w)] 
			(either (dim/1 * dim/2) = length? mdata: reduce mdata [
				m: make mtx [rows: dim/1 cols: dim/2 data: mdata]
			][
				cause-error 'user 'message ["Data length does not match dimensions!"]
			]
			)
		|	set w word! if (object? get/any w)(set w m: make mtx get w)
		| 	set m number!
		]	(insert matrices m)]
		unary-rule: [
			set unary [
				'transpose 
			| 	'rotate set n integer! 
			| 	'swap copy dims 3 skip 
			| 	'determinant
			|	'trace
			| 	'invert
			| 	'rref
			|	'identity opt [set d ['l | 'r]] 
			|	'game-of-life
			](
				insert unaries switch/default unary [
					rotate [reduce [unary n]] 
					swap [reduce [unary dims]]
					identity [either d [reduce [unary d]][unary]]
				][unary] 
			) expr-rule (
				unary: take unaries 
				switch unary [
					rotate [n: take unaries]
					swap [dims: take unaries]
					identity [if find ['l 'r] first unaries [d: take unaries]]
				]
				switch/default unary [
					rotate [self/rotate n matrices/1]
					;swap [matrices/1/(to-word rejoin ["swap-" dims/1]) dims/2 dims/3]
					swap [self/(to-word rejoin ["swap-" dims/1]) dims/2 dims/3 matrices/1]
					trace [insert matrices self/trace matrices/1]
					determinant [insert matrices self/determinant matrices/1]
					identity [insert matrices either d [
						self/identity/side matrices/1 d
					][
						self/identity matrices/1
					]]
					rref [insert matrices self/rref matrices/1]
					invert [insert matrices self/invert matrices/1]
				][
					self/:unary matrices/1
				]
			)
		]
		op-probe: [ahead [[pair! [block! | word!] | word!] ops-rule]]
		op-rule: [
			matrix-rule set op' ops-rule (insert ops op') expr-rule (
				op': take ops set [m2 m1] take/part matrices 2
				case [
					op' = 'div [op: :/ either number? m1 [m1: to-float m1][m1/to-float]]
					find [x augment] op' []
					true [op: get op']
				]
				case [
					all [number? reduce m1 number? reduce m2] [m1: (reduce m1) op reduce m2]
					number? reduce m1 [data: m2/data forall data [data/1: (reduce m1) op data/1] m1: m2]
					number? reduce m2 [data: m1/data forall data [data/1: data/1 op reduce m2]]
					true [case [
						find exclude ops-rule ['x 'augment] op' [m1: self/add :op m1 m2]
						(same? op' 'x) or (same? op' '×) [m1: self/multi m1 m2]
						same? op' 'X [m1: self/kronecker m1 m2]
						op' = 'augment [m1: self/augment m1 m2]
					]]
				]
				insert matrices m1
			)
		]
		expr-rule: [
			set var set-word! (insert vars var) 
			expr-rule (var: take vars set var copy matrices/1)
		|	ahead block! into rule 
		|	op-probe op-rule
		|	unary-rule
		| 	matrix-rule
		]
		parse spec rule: [some [
			ahead paren! into rule
		| 	expr-rule
		|	s: print ["No rule applied at: " :s]
		]]
		
		either number? m1: take matrices [m1][new-line/skip copy m1/data true m1/cols]
	]
]

Features:

• binary ops: ['+ | '- | '* | '/ | '% | '** | '>> | '<< | '>>> | 'and | 'or | 'xor | 'div | 'x]
• unary-matrix ops transpose and rotate n
• div turns args to floats
• x is standard matrix multiplication op
• * is Hadamard multiplication op
• X is Kronecker's multiplication op
• ops with scalar args
• order of ops as usual
• parens to change priority
• set-words to catch states

do %matrix.red
>> matrix [2x3 [1 2 3 4 5 6]]
== [
    1 2 3 
    4 5 6
]
;### Unary ops ###
; a) Transpose 
>> matrix [transpose 2x3 [1 2 3 4 5 6]]
== [
    1 4 
    2 5 
    3 6
]
; b) Rotate (clockwise [1 | 2 | 3]) 
>> matrix [rotate 1 2x3[1 2 3 4 5 6]]
== [
    4 1 
    5 2 
    6 3
]
>> matrix [rotate 2 2x3[1 2 3 4 5 6]]
== [
    6 5 4 
    3 2 1
]
>> matrix [rotate 3 2x3[1 2 3 4 5 6]]
== [
    3 6 
    2 5 
    1 4
]
;### Hadamard-like ops (preserving dims) ###
matrix [2x3 [1 2 3 4 5 6] + 2x3 [2 3 4 5 6 7]]
== [
    3 5 7 
    9 11 13
]
matrix [2x3 [1 2 3 4 5 6] * 2x3 [2 3 4 5 6 7]]
== [
    2 6 12 
    20 30 42
]
;### Standard product ###
matrix [2x3 [1 2 3 4 5 6] x 3x2 [2 3 4 5 6 7]]
== [
    28 34 
    64 79
]
;### Kronecker product ###
matrix [2x3 [1 2 3 4 5 6] X 3x2 [2 3 4 5 6 7]]
== [
    2 3 4 6 6 9 
    4 5 8 10 12 15 
    6 7 12 14 18 21 
    8 12 10 15 12 ...
;### Boolean logic ###
matrix [2x2 [1 0 1 1] and 2x2 [0 1 1 0]]
== [
    0 0 
    1 0
]
;### Scalar args ###
matrix [3x3 [1 2 3 2 4 5 3 6 2] * 2 - 5]
== [
    -3 -1 1 
    -1 3 5 
    1 7 -1
]
;### Prepared data ###
>> data1: [1 2 3 4 5 6]
== [1 2 3 4 5 6]
>> matrix [2x3 data1]
== [
    1 2 3 
    4 5 6
]
>> matrix [2x3 data1 x 3x3 [1 2 2 1 3 1 2 2 3]]
== [
    9 14 13 
    21 35 31
]
>> data2: object [rows: 3 cols: 2 data: [1 2 3 3 2 1]]
== make object! [
    rows: 3
    cols: 2
    data: [1 2 3 3 2 1]
]
>> matrix [2x3 data1 x data2]
== [
    13 11 
    31 29
]
;### Catching state ###
>> matrix [m: 3x3[1 2 3 4 5 6 7 8 9]]
== [
    1 2 3 
    4 5 6 
    7 8 9
]
>> probe m
make object! [
    rows: 3
    cols: 3
    data: [1 2 3 4 5 6 7 8 9]
    get-col: func [col][extract at data col cols]
    get-row: func [row][copy/part at data row - 1 * cols + 1 cols]
    get-idx: func [row col][pick data row - 1 * cols + col]
    to-float: func [][forall data [data/1: system/words/to-float data/1]]
    swap-dim: func [][c: cols cols: rows rows: c]
]
;### Pair entries ###
 matrix [2x2[1x1 1x2 2x1 2x2] x 2x2[1x2 2x1 3x2 2x3]]
== [
    4x6 4x7 
    8x6 8x7
]
;### Date entries ###
matrix [2x2[1-2-2017 1-3-2018 1-4-2017 31-5-2017] + 2x2[2 1 3 3]]
== [
    3-Feb-2017 2-Mar-2018 
    4-Apr-2017 3-Jun-2017
]
;### Different entries ###
>> matrix [2x2[1.2.3.4 1-3-2018 12:05 1%] + 2x2[2 1 60 3%]]
== [
    3.4.5.6 2-Mar-2018 
    12:06:00 4%
]
>> date: 2017-3-25 matrix [2x2[1.2.0.4 1-3-2018 12:05 1%] + 2x2[191.166.1.46 date/day 55 * 60 3%]]
== [
    192.168.1.50 26-Mar-2018 
    13:00:00 4%
]