2
2
MVEE := proc(XY,
3
3
{tolerance::positive:= 1e-4 }, #Convergence Criterion
4
4
{maxiterations::posint := 100 },
5
- {output::{ identical (data,plot), list ( identical (data,plot))} := data },
5
+ {showplot::truefalse := false },
6
6
{filled::truefalse := false }
7
7
)
8
8
9
- local alpha, evalues, evectors, i, l_error, ldata, ldataext, M, maxvalindex, n, ncols, nrows, p1, semiaxes, stepsize, U, U1, x, X, y;
10
- local A, center, l_output; #Output
11
-
12
- if hastype (output, ' list' ) then
13
- l_output := output;
14
- else
15
- l_output := [output];
16
- end if ;
9
+ local alpha, count, evalues, evectors, i, l_error, ldata, ldataext, M, maxvalindex, n, ncols, nrows, semiaxes, stepsize, U, U1, x, X, y;
10
+ local A, center; #Output
17
11
18
12
kernelopts (opaquemodules=false ):
19
13
20
14
ldata := Statistics:-PreProcessData(XY, 2 , ' copy' );
21
15
22
16
nrows, ncols := upperbound(ldata);
23
- ldataext := Matrix([ldata, Vector[column](nrows, ' :- fill' = 1 )], ' datatype = float' );
17
+ ldataext := Matrix([ldata, Vector[column](nrows, ' fill' = 1 )], ' datatype = float' );
24
18
25
19
if ncols <> 2 then
26
20
error " expected 2 columns of data, got %1" , ncols;
27
21
end if ;
28
22
23
+ count := 1 ;
29
24
l_error := 1 ;
30
25
31
26
U := Vector[column](1 ..nrows, ' fill' = 1 /nrows);
40
35
U1 := (1 - stepsize) * U;
41
36
U1[maxvalindex] := U1[maxvalindex] + stepsize;
42
37
l_error := LinearAlgebra:-Norm(LinearAlgebra:-DiagonalMatrix(U1 - U));
38
+ count := count + 1 ;
43
39
U := U1;
44
40
45
41
end do ;
46
42
47
43
A := (1 /ncols) * LinearAlgebra:-MatrixInverse(LinearAlgebra:-Transpose(ldata) . LinearAlgebra:-DiagonalMatrix(U) . ldata - (LinearAlgebra:-Transpose(ldata) . U) . LinearAlgebra:-Transpose((LinearAlgebra:-Transpose(ldata) . U)));
48
44
center := LinearAlgebra:-Transpose(ldata) . U;
49
45
evalues, evectors := LinearAlgebra:-Eigenvectors(A);
50
- evectors := evectors(.., sort[index](1 /~ (sqrt~(Re~(evalues))), `>`, ' :- output' = ' :- permutation' ));
46
+ evectors := evectors(.., sort[index](1 /~ (sqrt~(Re~(evalues))), `>`, ' output' = ' permutation' ));
51
47
semiaxes := sort(1 /~ (sqrt~(Re~(evalues))), `>`);
52
48
alpha := arctan(Re(evectors[2 ,1 ]) / Re(evectors[1 ,1 ]));
49
+ if showplot = true then
53
50
54
- if l_output = [' :-data' ] then
51
+ x := t -> center[1 ] + semiaxes[1 ] * cos (t) * cos(alpha) - semiaxes[2] * sin(t) * sin(alpha);
52
+ y := t -> center[2 ] + semiaxes[1 ] * cos (t) * sin(alpha) + semiaxes[2] * sin(t) * cos(alpha);
53
+ if filled then
54
+ return plots:-display(subs(CURVES=POLYGONS, plot([x (t), y (t), t = 0 ..2 *Pi], ' transparency' = 0.95 , _rest)));
55
+ else
56
+ return plot([x (t), y (t), t = 0 ..2 *Pi], _rest);
57
+
58
+ end if ;
59
+ else
55
60
return A, center, evectors, evalues;
56
- elif has ( l_output, ' :-plot' ) then
57
- x := t -> center[1] + semiaxes[1] * cos(t) * cos(alpha) - semiaxes[2] * sin(t) * sin(alpha);
58
- y := t -> center[2 ] + semiaxes[1 ] * cos (t) * sin(alpha) + semiaxes[2] * sin(t) * cos(alpha);
59
- if filled then
60
- p1 := plots:-display(subs(CURVES=POLYGONS, plot([x (t), y (t), t = 0 ..2 *Pi], ' :-transparency' = 0.95 , _rest)));
61
- else
62
- p1 := plot([x (t), y (t), t = 0 ..2 *Pi], _rest);
63
- end if ;
64
- return p1, `if `( has(l_output, ' :-data' ), op([A, center, evectors, evalues]), NULL );
65
61
end if ;
66
62
67
63
end proc:
0 commit comments