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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2671 - (download) (as text) (annotate)
Sun Aug 22 21:18:06 2004 UTC (15 years, 5 months ago) by dpvc
File size: 11028 byte(s)
Fixed some inconsistencies between handing of matrices within the
parser and Value packages.  Added a predefined Matrix context.

    1 ###########################################################################
    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' => \&Value::cmp,
   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->numberMatrix(@{$M}) if Value::isNumber($M->[0]);
   41   return $self->matrixMatrix(@{$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 $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     $x = Value::Real->make($x) if !ref($x);
   77     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 = (@dr,1); $r = $pkg->make($r)->transpose}
  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   return $self->TeX(undef,$self->{open},$self->{close}) if $$Value::context->flag('StringifyAsTeX');
  341   return $self->string(undef,$self->{open},$self->{close})
  342 }
  343 
  344 sub string {
  345   my $self = shift; my $equation = shift;
  346   my $def = ($equation->{context} || $$Value::context)->lists->get('Matrix');
  347   my $open = shift || $def->{open}; my $close = shift || $def->{close};
  348   my @coords = ();
  349   foreach my $x (@{$self->data}) {
  350     if (Value::isValue($x)) {push(@coords,$x->string($equation,$open,$close))}
  351       else {push(@coords,$x)}
  352   }
  353   return $open.join(',',@coords).$close;
  354 }
  355 
  356 #
  357 #  Use \matrix to lay out matrices
  358 #
  359 sub TeX {
  360   my $self = shift; my $equation = shift;
  361   my $def = ($equation->{context} || $$Value::context)->lists->get('Matrix');
  362   my $open = shift || $def->{open}; my $close = shift || $def->{close};
  363   $open = '\{' if $open eq '{'; $close = '\}' if $close eq '}';
  364   my $TeX = ''; my @entries = (); my $d;
  365   if ($self->isRow) {
  366     foreach my $x (@{$self->data}) {
  367       push(@entries,(Value::isValue($x))? $x->TeX($equation): $x);
  368     }
  369     $TeX .= join(' &',@entries) . "\n";
  370     $d = scalar(@entries);
  371   } else {
  372     foreach my $row (@{$self->data}) {
  373       foreach my $x (@{$row->data}) {
  374         push(@entries,(Value::isValue($x))? $x->TeX($equation): $x);
  375       }
  376       $TeX .= join(' &',@entries) . '\cr'."\n";
  377       $d = scalar(@entries); @entries = ();
  378     }
  379   }
  380   return '\left'.$open.'\begin{array}{'.('c'x$d).'}'."\n".$TeX.'\end{array}\right'.$close;
  381 }
  382 
  383 ###########################################################################
  384 
  385 1;
  386 

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9