Uh oh!
There was an error while loading. Please reload this page.
This repository was archived by the owner on Nov 20, 2020. It is now read-only.
forked from carvalho/numlua
- Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmatrix.lua
More file actions
Latest commit
374 lines (324 loc) · 9.87 KB
/
Copy pathmatrix.lua
File metadata and controls
374 lines (324 loc) · 9.87 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
--[[
-- matrix.lua
-- Multidimensional matrix library for NumericLua
-- Luis Carvalho (lexcarvalho@gmail.com)
-- See Copyright Notice in numlua.h
--]]
-- Methods
localnew, copy=matrix.new, matrix.copy
localsum, diag=matrix.sum, matrix.diag
localget, set, size=matrix.get, matrix.set, matrix.size
localconcat, shape=matrix.concat, matrix.shape
localtype, eps=numlua.type, mathx.eps
localunpack, assert, ipairs=unpack, assert, ipairs
localsetmetatable, select=setmetatable, select
localfloor, max=math.floor, math.max
localfunctioncheckmatrix (m)
localt=type(m)
assert(t=="matrix", "matrix expected, got " ..t)
returnm
end
localtranspose=matrix.transpose
matrix.t=transpose-- handy alias
functionmatrix.ctranspose (m) returntranspose(m, true) end
localzeros=function (...) returnset(new(...), 0) end
functionmatrix.ones (...) returnset(new(...), 1) end
functionmatrix.eye (n, c) returnset(set(new(n, n, c), 0), "D", 1) end
matrix.zeros=zeros
-- [ Metamethods ]
localmt=getmetatable(new(1))
localadd, mul, mmul=matrix.add, matrix.mul, matrix.mmul
localdiv, ls=matrix.div, matrix.ls
localsection, slice=matrix.section, matrix.slice
mt.__add=function (a, b)
iftype(a) =="number" ortype(a) =="complex" then
returnadd(b, a)
end
returnadd(a, b)
end
mt.__sub=function (a, b)
iftype(a) =="number" ortype(a) =="complex" then
returnadd(-b, a)
end
iftype(b) =="number" ortype(b) =="complex" then
returnadd(a, -b)
end
returnadd(a, b, -1)
end
mt.__mul=function (a, b)
iftype(a) =="number" ortype(a) =="complex" then
returnmul(b, a)
end
iftype(b) =="number" ortype(b) =="complex" then
returnmul(a, b)
end
localda, db=size(a, "#"), size(b, "#")
localn, m=size(a, 1), size(b, 2)
localiscomplex=a:iscomplex() orb:iscomplex()
ifda==1anddb==1then-- outer product?
returnmmul(zeros(n, n, iscomplex), a, b)
end
ifda==1then-- v * A?
returnmmul(zeros(m, iscomplex), b, a, "T")
end
ifdb==1then-- A * v?
returnmmul(zeros(n, iscomplex), a, b)
end
-- da = db = 2:
returnmmul(zeros(n, m, iscomplex), a, b)
end
mt.__mod=ls
mt.__div=function (a, b)
iftype(a) =="number" ortype(a) =="complex" then
returndiv(b, a, true)
end
iftype(b) =="number" ortype(b) =="complex" then
returndiv(a, b)
end
localx=ls(transpose(b), transpose(a))
returnsize(x, "#") ==2andtranspose(x) orx
end
-- TODO: __call using string triplets (based on section)
-- [ _section_(m, "f1:l1:s1, f2:l2:s2, ...") <=> m[[f1:l1:s1,...]] ]
-- ^----'triplet string'---^
mt.__call=function (a, ...)
localt=select(1, ...)
iftype(t) =="table" thenreturnsection(a, t) end
returnslice(a, ...)
end
localcabs, linspace=complex.abs, matrix.linspace
functionmatrix.seq (a, b, step)
locals=stepor1
localn=floor(cabs((b-a) /s+1))
returnlinspace(a, b, n)
end
functionmatrix.trace (m)
localr, c=shape(checkmatrix(m))
assert(size(m, "#") ==2andr==c, "square matrix expected")
returnsum(diag(m))
end
-- [ Logical ]
localfind, ifelse, which=matrix.find, matrix.ifelse, matrix.which
functionmatrix.any (m, cond) returnfind(m, cond) ~=nilend
functionmatrix.all (m, cond) returnfind(m, cond, true) ==nilend
-- count(m, cond) <=> fold(m, \a,e(a + (cond(e) and 1 or 0)), 0)
-- <=> sum(ifelse(copy(m), cond, 1, 0))
functionmatrix.count (m, cond) returnwhich(m, cond, "#") end
functionmatrix.merge (x, y, mask) returnifelse(copy(mask), 1, x, y) end
functionmatrix.pack (m, mask) returnwhich(m, mask, "v") end
functionmatrix.unpack (v, mask, m) returnset(m, which(m, mask), v) end
-- [ From/To table conversions ]
localfunctioncheckvector (t, iscomplex)
localisvector, iscomplex=true, iscomplexort.complex
fori, vinipairs(t) do
iftype(v) ~="number" andtype(v) ~="complex" then
isvector=false
iftype(v) =="matrix" andsize(v, "#") ==1then
v=v[1]
t[i] =v
isvector=true
end
end
ifisvectorthen
iscomplex=iscomplexortype(v) =="complex"
else
break
end
end
returnisvector, iscomplex
end
localfunctionfromtable (t, iscomplex)
assert(type(t) =="table", "table expected")
localisvector, iscomplex=checkvector(t, iscomplex)
ifisvectorthen-- base case?
localv=new(#t, iscomplex)
fori, einipairs(t) dov[i] =eend
returnv
end
-- recursion
fori, vinipairs(t) do
iftype(v) =="table" then-- recurse?
t[i] =fromtable(v, iscomplex)
end
end
-- fix if complex
iscomplex=false
for_, vinipairs(t) do
iscomplex=iscomplexorv:iscomplex()
end
ifiscomplexthen
fori, vinipairs(t) do
ifnotv:iscomplex() then
t[i] =v:complex()
end
end
end
returnconcat(unpack(t))
end
matrix.fromtable=fromtable
localfunctiontotable (m)
assert(type(m) =="matrix", "matrix expected")
locald, t=size(m, "#"), {}
fori=1, #mdo
t[i] =d==1andm[i] ortotable(m[i])
end
returnt
end
matrix.totable=totable
functionmatrix.list (m)
checkmatrix(m)
fori, einm:entries(true) do
localt= {m:eindex(i)}
t[#t+1] =e
print(unpack(t))
end
end
localfunctionformatnumber (x, d)
localfmt=dand ("%." ..d.."f") or"%g"
returnfmt:format(x)
end
localsignbit=mathx.signbit
localfunctionformatcomplex (c, d)
localre, im=c:real(), c:imag()
localfmt=signbit(im) and"%s%si" or"%s+%si"
returnfmt:format(formatnumber(re, d), formatnumber(im, d))
end
localfunctiongetmaxlen (fmt, d)
returnfunction (l, e) returnmax(l, #fmt(e, d)) end
end
localtconcat=table.concat
localfunctionprettyaux (v, ml, fmt, d) -- print vector with max length ml
localt= {}
fori=1, #vdo
localvi=fmt(v[i], d)
t[i] = (""):rep(3+ml-#vi) ..vi
end
returntconcat(t)
end
functionmatrix.pretty (m, d) -- `d` is number of decimal places
assert(size(checkmatrix(m), "#") <=2, "two-dimensional matrix expected")
localfmt=m:iscomplex() andformatcomplexorformatnumber
localml=m:fold(getmaxlen(fmt, d), 0) -- max length
ifsize(m, "#") ==1then
returnprettyaux(m, ml, fmt, d)
else-- m:size"#" == 2
localt= {}
fori=1, #mdot[i] =prettyaux(m[i], ml, fmt, d) end
returntconcat(t, "\n")
end
end
-- set metatable for class
matrix=setmetatable(matrix, {
__call=function(_, ...)
returntype(select(1, ...)) =="table" andfromtable(...) ornew(...)
end
})
-- [ Aggregators ]
localfunctionopfold (f, init)
localc
returnfunction (i, e)
ifi==1thenc=initend
c=f(c, v)
returnc
end
end
localsum2=function(x, y) returnx+yend
functionmatrix.cumsum (m)
returnm:apply(opfold(sum2, 0), true)
end
localprod2=function(x, y) returnx*yend
functionmatrix.cumprod (m)
returnm:apply(opfold(prod2, 1), true)
end
localprod=function (m) returnm:fold(prod2, 1) end
matrix.prod=prod
-- [ Linear algebra ]
localchol, lu, svd=matrix.chol, matrix.lu, matrix.svd
functionmatrix.kronecker (a, b)
assert(size(a, "#") ==2andsize(b, "#") ==2,
"two-dimensional matrix expected")
localra, ca, ica=shape(a, 1, true)
localrb, cb, icb=shape(b, 1, true)
localiscomplex=icaoricb
ifiscomplexthen
ifnoticathena=a:complex() end
ifnoticbthenb=b:complex() end
end
localx=new(ra*rb, ca*cb, iscomplex)
localindexr, indexc= {}, {}
localindex= {indexr, indexc}
fori=1, rado
localai=a[i]
indexr[1], indexr[2] = (i-1) *rb+1, i*rb
forj=1, cado
indexc[1], indexc[2] = (j-1) *cb+1, j*cb
mul(set(section(x, index), b), ai[j], true) -- x[index] = a[i][j] * b
end
end
returnx
end
functionmatrix.isposdef (m)
localc, msg=chol(checkmatrix(m))
ifc==nilthenerror(msg) end
returnnotc==false
end
functionmatrix.det (m)
localc=assert(lu(copy(checkmatrix(m)), true))
returnprod(diag(c))
end
functionmatrix.cond (m)
locals=assert(svd(checkmatrix(m), "n")) -- just singular values
returns[1] /s[#m]
end
-- effective rank from singular values `s`, max dim `m`, tolerance `tol`
locallt=function (x) returnfunction(e) returne<xendend
localfunctionsrank (s, m, tol)
localtol=tolor0
iftol<=0then-- set default tolerance?
tol=m*eps*s[1]
end
localr=s:find(lt(tol))
returnrandr-1or#s
end
functionmatrix.rank (m, tol)
locals=assert(svd(checkmatrix(m), "n")) -- just singular values
returnsrank(s, max(shape(m)), tol)
end
functionmatrix.null (m, tol)
localu, s, vh=assert(svd(checkmatrix(m)))
localnr, nc=shape(m)
localrank=srank(s, max(nr, nc), tol)
returnrank<ncandslice(vh, rank+1) ornil
end
functionmatrix.orth (m, tol)
localu=copy(checkmatrix(m))
locals=assert(svd(u, "l"))
localrank=srank(s, max(shape(m)), tol)
returnu{{}, {1, rank}} -- columns from 1 to rank
end
-- pseudo-inverse
functionmatrix.pinv (m, tol)
localu, s, vh=assert(svd(checkmatrix(m)))
localnr, nc=shape(m)
localrank=srank(s, max(nr, nc), tol)
localv=slice(vh, 1, rank)
fori=1, rankdo-- inv(s) * vh
v[i]:div(s[i], false, true) -- v[i,:] = v[i,:] / s[i], in-place
end
returnzeros(nc, nr):mmul(v, u{{}, {1, rank}}, "c", "c")
end
-- basic LS linear model fitting
functionmatrix.lm (a, b, svd)
localm, n=shape(checkmatrix(a))
assert(m>=n, "system is underdetermined")
assert(checkmatrix(b):size"#" ==1, "single RHS expected")
localx, rank=ls(a, b, svd)
-- report summary statistics
localcoef=slice(x, 1, n)
localrss= (b-a%coef):norm() ^2
localrss0= (b-b:sum() /m):norm() ^2
localdf=m-rank
localF=df/ (rank-1) * (rss0/rss-1)
localpvalue=1-stat.pf(F, rank-1, df)
return {coef=coef, rss=rss, df=df, F=F, pvalue=pvalue}
end