← Back to branch summary

~ubuntu-branches/debian/sid/genius/sid

« back to all changes in this revision

Viewing changes to lib/linear_algebra/misc.gel

  • Committer: Bazaar Package Importer
  • Author(s): Daniel Holbach
  • Date: 2006-08-21 12:57:45 UTC
  • Revision ID: james.westby@ubuntu.com-20060821125745-sl9ks8v7fq324bdf
Tags: upstream-0.7.6.1
Import upstream version 0.7.6.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/linear_algebra/misc.gel
 
1
SetHelp("ApplyOverMatrix", "matrix", "Apply a function over all entries of a matrix and return a matrix of the results")
 
2
function ApplyOverMatrix(a,func) = (
 
3
        if(not IsMatrix(a)) then
 
4
                (error("ApplyOverMatrix: argument 1 must be a matrix");bailout)
 
5
        else if(not IsFunction(func)) then
 
6
                (error("ApplyOverMatrix: argument 2 must be a function");bailout);
 
7
        for i = 1 to rows(a) do (
 
8
                for j = 1 to columns(a) do (
 
9
                        r@(i,j) = func(a@(i,j))
 
10
                )
 
11
        );
 
12
        r
 
13
);
 
14
protect("ApplyOverMatrix")
 
15
 
 
16
SetHelp("ApplyOverMatrix2", "matrix", "Apply a function over all entries of 2 matrices (or 1 value and 1 matrix) and return a matrix of the results")
 
17
function ApplyOverMatrix2(a,b,func) = (
 
18
        if(not IsMatrix(a) and not IsMatrix(b)) then
 
19
                (error("ApplyOverMatrix2: argument 1 or 2 must be a matrix");bailout)
 
20
        else if(not IsFunction(func)) then
 
21
                (error("ApplyOverMatrix2: argument 3 must be a function");bailout)
 
22
 
 
23
        else if(IsMatrix(a) and IsMatrix(b) and
 
24
                (rows(a)!=rows(b) or columns(a)!=columns(b))) then
 
25
                (error("ApplyOverMatrix2: cannot apply a function over two matrixes of different sizes");bailout);
 
26
        
 
27
        for i = 1 to rows(a) do (
 
28
                for j = 1 to columns(a) do (
 
29
                        if(IsMatrix(a)) then
 
30
                                aa = a@(i,j)
 
31
                        else
 
32
                                aa = a;
 
33
                        if(IsMatrix(b)) then
 
34
                                bb = b@(i,j)
 
35
                        else
 
36
                                bb = b;
 
37
                        r@(i,j) = func(aa,bb)
 
38
                )
 
39
        );
 
40
        r
 
41
);
 
42
protect("ApplyOverMatrix2")
 
43
 
 
44
#calculate a trace function
 
45
SetHelp("Trace", "linear_algebra","Calculate the trace of a matrix");
 
46
function Trace(m) = (
 
47
        if(not IsMatrix(m) or not IsValueOnly(m)) then
 
48
                (error("Trace: argument not a value only matrix");bailout)
 
49
        else if(rows(m)!=columns(m)) then
 
50
                (error("Trace: matrix not a square");bailout);
 
51
        sum i = 1 to rows(m) do m@(i,i)
 
52
);
 
53
protect("Trace")
 
54
SetHelpAlias("Trace", "trace")
 
55
trace = Trace
 
56
protect("trace")
 
57
 
 
58
#calculate convolution of two horizontal vectors
 
59
SetHelp("Convolution", "linear_algebra","Calculate convolution of two horizontal vectors");
 
60
function Convolution(a,b) = (
 
61
        if(not IsMatrix(a) or not IsValueOnly(a) or
 
62
           not IsMatrix(b) or not IsValueOnly(b) or
 
63
           rows(a)>1 or rows(b)>1) then
 
64
                (error("Convolution: arguments not value only horizontal vectors");bailout)
 
65
        else if(columns(a)!=columns(b)) then
 
66
                (error("Convolution: arguments must be identical vectors");bailout);
 
67
        ca = columns(a);
 
68
        sum i = 1 to ca do a@(1,i)*b@(1,ca-i+1)
 
69
);
 
70
protect("Convolution")
 
71
SetHelpAlias("Convolution", "convol")
 
72
convol = Convolution
 
73
protect("convol")
 
74
 
 
75
#calculate convolution of two horizontal vectors and return the result
 
76
#not added together but in a vector
 
77
SetHelp("ConvolutionVector", "linear_algebra","Calculate convolution of two horizontal vectors");
 
78
function ConvolutionVector(a,b) = (
 
79
        if(not IsMatrix(a) or not IsValueOnly(a) or
 
80
           not IsMatrix(b) or not IsValueOnly(b) or
 
81
           rows(a)>1 or rows(b)>1) then
 
82
                (error("ConvolutionVector: arguments not value only horizontal vectors");bailout)
 
83
        else if(columns(a)!=columns(b)) then
 
84
                (error("ConvolutionVector: arguments must be identical vectors");bailout);
 
85
        ca = columns(a);
 
86
        r = zeros (1,ca);
 
87
        for i = 1 to ca do (
 
88
                r@(1,i) = a@(1,i)*b@(1,ca-i+1)
 
89
        );
 
90
        r
 
91
);
 
92
protect("ConvolutionVector")
 
93
 
 
94
#calculate the sum of all elements in a matrix
 
95
SetHelp("MatrixSum", "matrix","Calculate the sum of all elements in a matrix");
 
96
function MatrixSum(a) = (
 
97
        if(not IsMatrix(a) or not IsValueOnly(a)) then
 
98
                (error("MatrixSum: argument not a value only matrix");bailout);
 
99
        sum n in a do n
 
100
);
 
101
protect("MatrixSum")
 
102
 
 
103
SetHelp("MatrixSumSquares", "matrix","Calculate the sum of squares of all elements in a matrix");
 
104
function MatrixSumSquares(a) = (
 
105
        if(not IsMatrix(a) or not IsValueOnly(a)) then
 
106
                (error("MatrixSumSquares: argument not a value only matrix");bailout);
 
107
        sum n in a do n^2
 
108
);
 
109
protect("MatrixSumSquares")
 
110
 
 
111
#calculate the product of all elements in a matrix
 
112
SetHelp("MatrixProduct","matrix", "Calculate the product of all elements in a matrix")
 
113
function MatrixProduct(a) = (
 
114
        if(not IsMatrix(a) or not IsValueOnly(a)) then
 
115
                (error("matprod: argument not a value only matrix");bailout);
 
116
        prod n in a do n
 
117
);
 
118
protect("MatrixProduct")
 
119
 
 
120
SetHelp("Submatrix", "matrix", "Return column(s) and row(s) from a matrix")
 
121
function Submatrix(m,r,c) = [m@(r,c)]
 
122
protect ("Submatrix")
 
123
 
 
124
SetHelp("ComplementSubmatrix", "matrix", "Remove column(s) and row(s) from a matrix");
 
125
function ComplementSubmatrix(m,r,c) = [m@(IndexComplement(r, rows(m)), IndexComplement (c, columns (m)))]
 
126
protect ("ComplementSubmatrix")
 
127
 
 
128
# Minor of a matrix (determinant of a submatrix given by deleting
 
129
# one row and one column)
 
130
SetHelp("Minor","linear_algebra", "Get the i-j minor of a matrix")
 
131
function Minor(M,i,j) = det (ComplementSubmatrix (M, i, j))
 
132
protect("Minor");
 
133
 
 
134
#classical adjoint (adjugate) of a matrix
 
135
SetHelp("adj","linear_algebra", "Get the classical adjoint (adjugate) of a matrix");
 
136
function adj(m) = (
 
137
        if(not IsMatrix(m) or not IsValueOnly(m)) then
 
138
                (error("adj: argument not a value-only matrix");bailout)
 
139
        else if(rows(m)!=columns(m)) then
 
140
                (error("adj: argument not a square matrix");bailout)
 
141
        else if(rows(m)<2) then
 
142
                (error("adj: argument cannot be 1x1 matrix");bailout);
 
143
 
 
144
        a = zeros (rows(m),rows(m));
 
145
        for i = 1 to rows(m) do (
 
146
                for j = 1 to rows(m) do (
 
147
                        a@(j,i) = ((-1)^(i+j))*Minor(m,i,j);
 
148
                );
 
149
        );
 
150
        a
 
151
);
 
152
protect("adj")
 
153
SetHelpAlias ("adj", "Adjugate")
 
154
Adjugate = adj
 
155
protect("Adjugate")
 
156
 
 
157
SetHelp("MinimizeFunction","functions","Find the first value where f(x)=0");
 
158
function MinimizeFunction(func,x,incr) = (
 
159
        if(not IsValue(x) or not IsValue(incr)) then
 
160
                (error("MinimizeFunction: x,incr arguments not values");bailout)
 
161
        else if(not IsFunction(func)) then
 
162
                (error("MinimizeFunction: func argument not a function");bailout);
 
163
        while(func(x)>0) do x=x+incr;
 
164
        x
 
165
);
 
166
protect("MinimizeFunction")
 
167
 
 
168
SetHelp("MakeDiagonal","matrix","Make diagonal matrix from a vector");
 
169
function MakeDiagonal(v,arg...) = (
 
170
        if IsValue (v) and IsNull (arg) then
 
171
                return [v]
 
172
        else if IsMatrix (v) and IsNull (arg) then
 
173
                m = v
 
174
        else if IsValue (v) and IsMatrix (arg) then
 
175
                m = [v,arg]
 
176
        else
 
177
                (error("MakeDiagonal: arguments not a vector or a list of values");bailout);
 
178
        r = zeros (elements(m),elements(m));
 
179
        for i = 1 to elements(m) do (
 
180
                r@(i,i) = m@(i)
 
181
        );
 
182
        r
 
183
);
 
184
protect("MakeDiagonal")
 
185
SetHelpAlias("MakeDiagonal","diag")
 
186
diag = MakeDiagonal
 
187
protect("diag")
 
188
 
 
189
SetHelp("SwapRows","matrix","Swap two rows in a matrix");
 
190
function SwapRows(m,row1,row2) = (
 
191
        if(not IsMatrix(m) or not IsInteger(row1) or
 
192
           not IsInteger(row2)) then
 
193
                (error("SwapRows: arguments are not the right type");bailout)
 
194
        else if(row1>rows(m) or row2>rows(m)) then
 
195
                (error("SwapRows: argument out of range");bailout)
 
196
        else if(row1 != row2) then (
 
197
                tmp = m@(row1,);
 
198
                m@(row1,) = m@(row2,);
 
199
                m@(row2,) = tmp
 
200
        );
 
201
        m
 
202
);
 
203
protect("SwapRows")
 
204
 
 
205
SetHelp("RowSum","matrix","Calculate sum of each row in a matrix");
 
206
function RowSum(m) = (
 
207
        if not IsMatrix(m) then
 
208
                (error("RowSum: argument not matrix");bailout);
 
209
        r = zeros (rows(m),1);
 
210
        for i = 1 to rows(m) do (
 
211
                for j = 1 to columns(m) do
 
212
                        r@(i,1) = r@(i,1) + m@(i,j)
 
213
        );
 
214
        r
 
215
);
 
216
protect("RowSum")
 
217
 
 
218
SetHelp("RowSumSquares","matrix","Calculate sum of squares of each row in a matrix");
 
219
function RowSumSquares(m) = (
 
220
        if not IsMatrix(m) then
 
221
                (error("RowSumSquares: argument not matrix");bailout);
 
222
        r = zeros (rows(m),1);
 
223
        for i = 1 to rows(m) do (
 
224
                for j = 1 to columns(m) do
 
225
                        r@(i,1) = r@(i,1) + m@(i,j)^2
 
226
        );
 
227
        r
 
228
);
 
229
protect("RowSumSquares")
 
230
 
 
231
#sort a horizontal vector
 
232
SetHelp("SortVector","matrix","Sort vector elements");
 
233
function SortVector(v) = (
 
234
        if not IsVector(v) then
 
235
                (error("SortVector: argument not a vector");bailout);
 
236
        j = 1;
 
237
        do (
 
238
                sorted = true;
 
239
                for i = 1 to elements(v)-j do (
 
240
                        if v@(i)>v@(i+1) then (
 
241
                                t = v@(i);
 
242
                                v@(i) = v@(i+1);
 
243
                                v@(i+1) = t;
 
244
                                sorted = false
 
245
                        )
 
246
                );
 
247
                j=j+1
 
248
        ) while not sorted;
 
249
        v
 
250
);
 
251
protect("SortVector")
 
252
 
 
253
SetHelp("ReverseVector","matrix","Reverse elements in a vector");
 
254
function ReverseVector(v) = (
 
255
        if not IsVector(v) then
 
256
                (error("ReverseVector: argument not a vector");bailout);
 
257
        r = zeros (rows(v),columns(v));
 
258
        ev = elements(v);
 
259
        for i = 1 to ev do
 
260
                r@(i) = v@(ev-i+1);
 
261
        r
 
262
);
 
263
protect("ReverseVector")
 
264
 
 
265
SetHelp("IsVector", "matrix", "Is argument a horizontal or a vertical vector")
 
266
function IsVector(v) = (IsMatrix(v) and (columns(v) == 1 or rows(v) == 1))
 
267
protect("IsVector")
 
268
 
 
269
SetHelp("UpperTriangular", "matrix", "Zero out entries below the diagonal")
 
270
function UpperTriangular(M) = (
 
271
        if not IsMatrix(M) or not IsMatrixSquare(M) then
 
272
                (error("UpperTriangular: argument not a square matrix");bailout);
 
273
        for i=2 to rows(M) do (
 
274
                for j=1 to i-1 do (
 
275
                        M@(i,j) = 0
 
276
                )
 
277
        );
 
278
        M
 
279
)
 
280
protect ("UpperTriangular")
 
281
 
 
282
SetHelp("LowerTriangular", "matrix", "Zero out entries above the diagonal")
 
283
function LowerTriangular(M) = (
 
284
        if not IsMatrix(M) or not IsMatrixSquare(M) then
 
285
                (error("LowerTriangular: argument not a square matrix");bailout);
 
286
        UpperTriangular (M.').'
 
287
)
 
288
protect ("LowerTriangular")
 
289
 
 
290
# FIXME: stupid
 
291
SetHelp("IsUpperTriangular", "matrix", "Is a matrix upper triangular")
 
292
function IsUpperTriangular(M) = (
 
293
        if not IsMatrix(M) or not IsMatrixSquare(M) then
 
294
                (error("IsUpperTriangular: argument not a square matrix");bailout);
 
295
        M == UpperTriangular(M)
 
296
)
 
297
protect ("IsUpperTriangular")
 
298
 
 
299
# FIXME: stupid
 
300
SetHelp("IsLowerTriangular", "matrix", "Is a matrix lower triangular")
 
301
function IsLowerTriangular(M) = (
 
302
        if not IsMatrix(M) or not IsMatrixSquare(M) then
 
303
                (error("IsLowerTriangular: argument not a square matrix");bailout);
 
304
        M == LowerTriangular(M)
 
305
)
 
306
protect ("IsLowerTriangular")
 
307
 
 
308
# FIXME: stupid
 
309
SetHelp("IsDiagonal", "matrix", "Is a matrix diagonal")
 
310
function IsDiagonal(M) = (
 
311
        if not IsMatrix(M) or not IsMatrixSquare(M) then
 
312
                (error("IsDiagonal: argument not a square matrix");bailout);
 
313
        M == UpperTriangular(LowerTriangular(M))
 
314
)
 
315
protect ("IsDiagonal")
 
316
 
 
317
SetHelp("CompoundMatrix", "matrix", "Calculate the kth compund matrix of A")
 
318
function CompoundMatrix(k,A) = (
 
319
        if not IsInteger(k) or k < 1 or k > min(columns(A),rows(A)) or not IsMatrix(A) then
 
320
                (error("CompoundMatrix: arguments of right type/size");bailout);
 
321
        C=[0];
 
322
        gamma = Combinations(k,rows(A));
 
323
        omega = Combinations(k,columns(A));
 
324
        for i=1 to elements(gamma) do
 
325
                for j=1 to elements(omega) do
 
326
                        C@(i,j) = det (A@(gamma@(i),omega@(j)));
 
327
        C
 
328
)
 
329
protect ("CompoundMatrix")
 
330
 
 
331
SetHelp("MakeVector", "matrix", "Make column vector out of matrix by putting columns above each other")
 
332
function MakeVector(A) = (
 
333
        if IsNull(A) then return null
 
334
        else if not IsMatrix(A) then
 
335
                (error("MakeVector: arguments not a matrix");bailout);
 
336
        r = null;
 
337
        for k=1 to columns(A) do (
 
338
                r = [r;A@(,k)]
 
339
        );
 
340
        r
 
341
)
 
342
protect ("MakeVector")