[system] / trunk / pg / lib / Value / Matrix.pm Repository:
ViewVC logotype

Annotation of /trunk/pg/lib/Value/Matrix.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2579 - (view) (download) (as text)

1 : sh002i 2558 ###########################################################################
2 :     #
3 :     # Implements the Matrix class.
4 :     #
5 :     # @@@ Still needs lots of work @@@
6 :     #
7 :     package Value::Matrix;
8 :     my $pkg = 'Value::Matrix';
9 :    
10 :     use strict;
11 :     use vars qw(@ISA);
12 :     @ISA = qw(Value);
13 :    
14 :     use overload
15 :     '+' => \&add,
16 :     '-' => \&sub,
17 :     '*' => \&mult,
18 :     '/' => \&div,
19 :     '**' => \&power,
20 :     '.' => \&Value::_dot,
21 :     'x' => \&Value::cross,
22 :     '<=>' => \&compare,
23 :     'cmp' => \&compare,
24 :     'neg' => sub {$_[0]->neg},
25 :     'nomethod' => \&Value::nomethod,
26 :     '""' => \&stringify;
27 :    
28 :     #
29 :     # Convert a value to a matrix. The value can be:
30 :     # a list of numbers or list of (nested) references to arrays of numbers
31 :     # a point, vector or matrix object
32 :     #
33 :     sub new {
34 :     my $self = shift; my $class = ref($self) || $self;
35 :     my $M = shift;
36 :     return bless {data => $M->data}, $class
37 :     if (Value::class($M) =~ m/Point|Vector|Matrix/ && scalar(@_) == 0);
38 :     $M = [$M,@_] if ((defined($M) && ref($M) ne 'ARRAY') || scalar(@_) > 0);
39 :     Value::Error("Matrices must have at least one entry") unless defined($M) && scalar(@{$M}) > 0;
40 :     return $self->matrixMatrix(@{$M}) if (ref($M->[0]) eq 'ARRAY');
41 :     return $self->numberMatrix(@{$M});
42 :     }
43 :    
44 :     #
45 :     # (Recusrively) make a matrix from a list of array refs
46 :     # and report errors about the entry types
47 :     #
48 :     sub matrixMatrix {
49 :     my $self = shift; my $class = ref($self) || $self;
50 :     my ($x,$m); my @M = (); my $d = scalar(@{$_[0]}); my $isFormula = 0;
51 :     foreach $x (@_) {
52 :     if (Value::isFormula($x)) {push(@M,$x); $isFormula = 1} else {
53 :     $m = $pkg->new($x); push(@M,$m);
54 :     $isFormula = 1 if Value::isFormula($m);
55 :     }
56 :     }
57 :     my ($type,$len) = ($M[0]->entryType->{name},$M[0]->length);
58 :     foreach $x (@M) {
59 :     Value::Error("Matrix rows must all be the same type")
60 :     unless (defined($x->entryType) && $type eq $x->entryType->{name});
61 :     Value::Error("Matrix rows must all be the same length") unless ($len eq $x->length);
62 :     }
63 :     return $self->formula([@M]) if $isFormula;
64 :     bless {data => [@M]}, $class;
65 :     }
66 :    
67 :     #
68 :     # Form a 1 x n matrix from a list of numbers
69 :     # (could become a row of an m x n matrix)
70 :     #
71 :     sub numberMatrix {
72 :     my $self = shift; my $class = ref($self) || $self;
73 :     my @M = (); my $isFormula = 0;
74 :     foreach my $x (@_) {
75 :     Value::Error("Matrix row entries must be numbers") unless (Value::isNumber($x));
76 : dpvc 2579 $x = Value::Real->make($x) if $$Value::context->flag('useFuzzyReals') && !Value::isFormula($x);
77 : sh002i 2558 push(@M,$x); $isFormula = 1 if Value::isFormula($x);
78 :     }
79 :     return $self->formula([@M]) if $isFormula;
80 :     bless {data => [@M]}, $class;
81 :     }
82 :    
83 :     #
84 :     # Recursively get the entries in the matrix and return
85 :     # an array of (references to arrays of ... ) numbers
86 :     #
87 :     sub value {
88 :     my $self = shift;
89 :     my $M = $self->data;
90 :     return @{$M} if Value::class($M->[0]) ne 'Matrix';
91 :     my @M = ();
92 :     foreach my $x (@{$M}) {push(@M,[$x->value])}
93 :     return @M;
94 :     }
95 :     #
96 :     # The number of rows in the matrix (for n x m)
97 :     # or the number of entries in a 1 x n matrix
98 :     #
99 :     sub length {return scalar(@{shift->{data}})}
100 :     #
101 :     # Recursively get the dimensions of the matrix.
102 :     # Returns (n) for a 1 x n, or (n,m) for an n x m, etc.
103 :     #
104 :     sub dimensions {
105 :     my $self = shift;
106 :     my $r = $self->length;
107 :     my $v = $self->data;
108 :     return ($r,) if (Value::class($v->[0]) ne 'Matrix');
109 :     return ($r,$v->[0]->dimensions);
110 :     }
111 :     #
112 :     # Return the proper type for the matrix
113 :     #
114 :     sub typeRef {
115 :     my $self = shift;
116 :     return Value::Type($self->class, $self->length, $Value::Type{number})
117 :     if (Value::class($self->data->[0]) ne 'Matrix');
118 :     return Value::Type($self->class, $self->length, $self->data->[0]->typeRef);
119 :     }
120 :    
121 :     #
122 :     # True if the matrix is a square matrix
123 :     #
124 :     sub isSquare {
125 :     my $self = shift;
126 :     my @d = $self->dimensions;
127 :     return 0 if scalar(@d) > 2;
128 :     return 1 if scalar(@d) == 1 && $d[0] == 1;
129 :     return $d[0] == $d[1];
130 :     }
131 :    
132 :     #
133 :     # True if the matrix is 1-dimensional (i.e., is a matrix row)
134 :     #
135 :     sub isRow {
136 :     my $self = shift;
137 :     my @d = $self->dimensions;
138 :     return scalar(@d) == 1;
139 :     }
140 :    
141 :     #
142 :     # Make arbitrary data into a matrix, if possible
143 :     #
144 :     sub promote {
145 :     my $x = shift;
146 :     return $pkg->new($x,@_) if scalar(@_) > 0 || ref($x) eq 'ARRAY';
147 :     return $x if ref($x) eq $pkg;
148 :     return $pkg->make(@{$x->data}) if Value::class($x) =~ m/Point|Vector/;
149 :     Value::Error("Can't convert ".Value::showClass($x)." to a Matrix");
150 :     }
151 :    
152 :     ############################################
153 :     #
154 :     # Operations on matrices
155 :     #
156 :    
157 :     sub add {
158 :     my ($l,$r,$flag) = @_;
159 :     if ($l->promotePrecedence($r)) {return $r->add($l,!$flag)}
160 :     ($l,$r) = (promote($l)->data,promote($r)->data);
161 :     Value::Error("Matrix addition with different dimensions")
162 :     unless scalar(@{$l}) == scalar(@{$r});
163 :     my @s = ();
164 :     foreach my $i (0..scalar(@{$l})-1) {push(@s,$l->[$i] + $r->[$i])}
165 :     return $pkg->make(@s);
166 :     }
167 :    
168 :     sub sub {
169 :     my ($l,$r,$flag) = @_;
170 :     if ($l->promotePrecedence($r)) {return $r->sub($l,!$flag)}
171 :     ($l,$r) = (promote($l)->data,promote($r)->data);
172 :     Value::Error("Matrix subtraction with different dimensions")
173 :     unless scalar(@{$l}) == scalar(@{$r});
174 :     if ($flag) {my $tmp = $l; $l = $r; $r = $tmp};
175 :     my @s = ();
176 :     foreach my $i (0..scalar(@{$l})-1) {push(@s,$l->[$i] - $r->[$i])}
177 :     return $pkg->make(@s);
178 :     }
179 :    
180 :     sub mult {
181 :     my ($l,$r,$flag) = @_;
182 :     if ($l->promotePrecedence($r)) {return $r->mult($l,!$flag)}
183 :     #
184 :     # Constant multiplication
185 :     #
186 :     if (Value::matchNumber($r) || Value::isComplex($r)) {
187 :     my @coords = ();
188 :     foreach my $x (@{$l->data}) {push(@coords,$x*$r)}
189 :     return $pkg->make(@coords);
190 :     }
191 :     #
192 :     # Make points and vectors into columns if they are on the right
193 :     #
194 :     if (!$flag && Value::class($r) =~ m/Point|Vector/)
195 :     {$r = (promote($r))->transpose} else {$r = promote($r)}
196 :     #
197 :     if ($flag) {my $tmp = $l; $l = $r; $r = $tmp}
198 :     my @dl = $l->dimensions; my @dr = $r->dimensions;
199 :     if (scalar(@dl) == 1) {@dl = (1,@dl); $l = $pkg->make($l)}
200 :     if (scalar(@dr) == 1) {@dr = (1,@dr); $r = $pkg->make($r)}
201 :     Value::Error("Can only multiply 2-dimensional matrices") if scalar(@dl) > 2 || scalar(@dr) > 2;
202 :     Value::Error("Matices of dimensions $dl[0]x$dl[1] and $dr[0]x$dr[1] can't be multiplied")
203 :     unless ($dl[1] == $dr[0]);
204 :     #
205 :     # Do matrix multiplication
206 :     #
207 :     my @l = $l->value; my @r = $r->value;
208 :     my @M = ();
209 :     foreach my $j (0..$dr[1]-1) {
210 :     my @row = ();
211 :     foreach my $i (0..$dl[0]-1) {
212 :     my $s = 0;
213 :     foreach my $k (0..$dl[1]-1) {$s += $l[$i]->[$k] * $r[$k]->[$j]}
214 :     push(@row,$s);
215 :     }
216 :     push(@M,$pkg->make(@row));
217 :     }
218 :     return $pkg->make(@M);
219 :     }
220 :    
221 :     sub div {
222 :     my ($l,$r,$flag) = @_;
223 :     if ($l->promotePrecedence($r)) {return $r->div($l,!$flag)}
224 :     Value::Error("Can't divide by a Matrix") if $flag;
225 :     Value::Error("Matrices can only be divided by numbers")
226 :     unless (Value::matchNumber($r) || Value::isComplex($r));
227 :     Value::Error("Division by zero") if $r == 0;
228 :     my @coords = ();
229 :     foreach my $x (@{$l->data}) {push(@coords,$x/$r)}
230 :     return $pkg->make(@coords);
231 :     }
232 :    
233 :     sub power {
234 :     my ($l,$r,$flag) = @_;
235 :     if ($l->promotePrecedence($r)) {return $r->power($l,!$flag)}
236 :     Value::Error("Can't use Matrices in exponents") if $flag;
237 :     Value::Error("Only square matrices can be raised to a power") unless $l->isSquare;
238 :     return Value::Matrix::I($l->length) if $r == 0;
239 :     Value::Error("Matrix powers must be positive integers") unless $r =~ m/^[1-9]\d*$/;
240 :     my $M = $l; foreach my $i (2..$r) {$M = $M*$l}
241 :     return $M;
242 :     }
243 :    
244 :     #
245 :     # Do lexicographic comparison
246 :     #
247 :     sub compare {
248 :     my ($l,$r,$flag) = @_;
249 :     if ($l->promotePrecedence($r)) {return $r->compare($l,!$flag)}
250 :     ($l,$r) = (promote($l)->data,promote($r)->data);
251 :     Value::Error("Matrix comparison with different dimensions")
252 :     unless scalar(@{$l}) == scalar(@{$r});
253 :     if ($flag) {my $tmp = $l; $l = $r; $r = $tmp};
254 :     my $cmp = 0;
255 :     foreach my $i (0..scalar(@{$l})-1) {
256 :     $cmp = $l->[$i] <=> $r->[$i];
257 :     last if $cmp;
258 :     }
259 :     return $cmp;
260 :     }
261 :    
262 :     sub neg {
263 :     my $p = promote(@_)->data;
264 :     my @coords = ();
265 :     foreach my $x (@{$p}) {push(@coords,-$x)}
266 :     return $pkg->make(@coords);
267 :     }
268 :    
269 :     #
270 :     # Transpose an n x m matrix
271 :     #
272 :     sub transpose {
273 :     my $self = shift;
274 :     my @d = $self->dimensions;
275 :     if (scalar(@d) == 1) {@d = (1,@d); $self = $pkg->make($self)}
276 :     Value::Error("Can't transpose ".scalar(@d)."-dimensional matrices") unless scalar(@d) == 2;
277 :     my @M = (); my $M = $self->data;
278 :     foreach my $j (0..$d[1]-1) {
279 :     my @row = ();
280 :     foreach my $i (0..$d[0]-1) {push(@row,$M->[$i]->data->[$j])}
281 :     push(@M,$pkg->make(@row));
282 :     }
283 :     return $pkg->make(@M);
284 :     }
285 :    
286 :     #
287 :     # Get an identity matrix of the requested size
288 :     #
289 :     sub I {
290 :     my $d = shift; $d = shift if ref($d) eq $pkg;
291 :     my @M = (); my @Z = split('',0 x $d);
292 :     foreach my $i (0..$d-1) {
293 :     my @row = @Z; $row[$i] = 1;
294 :     push(@M,$pkg->make(@row));
295 :     }
296 :     return $pkg->make(@M);
297 :     }
298 :    
299 :     #
300 :     # Extract a given row from the matrix
301 :     #
302 :     sub row {
303 :     my $M = promote(shift); my $i = shift;
304 :     return if $i == 0; $i-- if $i > 0;
305 :     if ($M->isRow) {return if $i != 0; return $M}
306 :     return $M->data->[$i];
307 :     }
308 :    
309 :     #
310 :     # Extract a given element from the matrix
311 :     #
312 :     sub element {
313 :     my $M = promote(shift);
314 :     return $M->extract(@_);
315 :     }
316 :    
317 :     #
318 :     # Extract a given column from the matrix
319 :     #
320 :     sub column {
321 :     my $M = promote(shift); my $j = shift;
322 :     return if $j == 0; $j-- if $j > 0;
323 :     my @d = $M->dimensions; my @col = ();
324 :     return if $j+1 > $d[1];
325 :     return $M->data->[$j] if scalar(@d) == 1;
326 :     foreach my $row (@{$M->data}) {push(@col,$pkg->make($row->data->[$j]))}
327 :     return $pkg->make(@col);
328 :     }
329 :    
330 :     # @@@ removeRow, removeColumn @@@
331 :     # @@@ Det, inverse @@@
332 :    
333 :     ############################################
334 :     #
335 :     # Generate the various output formats
336 :     #
337 :    
338 :     sub stringify {
339 :     my $self = shift;
340 : dpvc 2579 my $open = $$Value::context->lists->get('Matrix')->{open};
341 :     my $close = $$Value::context->lists->get('Matrix')->{close};
342 : sh002i 2558 return $open.join(',',@{$self->data}).$close
343 :     if (Value::class($self->data->[0]) ne 'Matrix');
344 :     return $open.join(",\n ",@{$self->data}).$close;
345 :     }
346 :    
347 :     sub string {
348 :     my $self = shift; my $equation = shift;
349 : dpvc 2579 my $open = shift || $$Value::context->lists->get('Matrix')->{open};
350 :     my $close = shift || $$Value::context->lists->get('Matrix')->{close};
351 : sh002i 2558 my @coords = ();
352 :     foreach my $x (@{$self->data}) {
353 :     if (Value::isValue($x)) {push(@coords,$x->string($equation,$open,$close))}
354 :     else {push(@coords,$x)}
355 :     }
356 :     return $open.join(',',@coords).$close;
357 :     }
358 :    
359 :     #
360 :     # Use \matrix to lay out matrices
361 :     #
362 :     sub TeX {
363 :     my $self = shift; my $equation = shift;
364 : dpvc 2579 my $open = shift || $$Value::context->lists->get('Matrix')->{open};
365 :     my $close = shift || $$Value::context->lists->get('Matrix')->{close};
366 : sh002i 2558 $open = '\{' if $open eq '{'; $close = '\}' if $close eq '}';
367 :     my $TeX = ''; my @entries = ();
368 : dpvc 2579 if ($self->isRow) {
369 :     foreach my $x (@{$self->data}) {
370 : sh002i 2558 push(@entries,(Value::isValue($x))? $x->TeX($equation,$open,$close): $x);
371 :     }
372 : dpvc 2579 $TeX .= join(' &',@entries) . "\n";
373 :     } else {
374 :     foreach my $row (@{$self->data}) {
375 :     foreach my $x (@{$row->data}) {
376 :     push(@entries,(Value::isValue($x))? $x->TeX($equation,$open,$close): $x);
377 :     }
378 :     $TeX .= join(' &',@entries) . '\cr'."\n"; @entries = ();
379 :     }
380 : sh002i 2558 }
381 :     return '\left'.$open.'\matrix{'."\n".$TeX.'}\right'.$close;
382 :     }
383 :    
384 :     ###########################################################################
385 :    
386 :     1;
387 :    

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9