[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 2800 - (download) (as text) (annotate)
Sun Sep 19 14:27:39 2004 UTC (15 years, 5 months ago) by dpvc
File size: 11496 byte(s)
Added isZero and isOne checks for Parser::Value objects (i.e., for
constants within formulas).  These now correctly handle vector and
matrices, in particular.  The isOne and isZero checks are used in the
reduce() method to simplify formulas.

    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 #  See if the matrix is an Indenity matrix
  143 #
  144 sub isOne {
  145   my $self = shift;
  146   return 0 unless $self->isSquare;
  147   my $i = 0;
  148   foreach my $row (@{$self->{data}}) {
  149     my $j = 0;
  150     foreach my $k (@{$row->{data}}) {
  151       return 0 unless $k eq (($i == $j)? "1": "0");
  152       $j++;
  153     }
  154     $i++;
  155   }
  156   return 1;
  157 }
  158 
  159 #
  160 #  See if the matrix is all zeros
  161 #
  162 sub isZero {
  163   my $self = shift;
  164   foreach my $x (@{$self->{data}}) {return 0 unless $x->isZero}
  165   return 1;
  166 }
  167 
  168 #
  169 #  Make arbitrary data into a matrix, if possible
  170 #
  171 sub promote {
  172   my $x = shift;
  173   return $pkg->new($x,@_) if scalar(@_) > 0 || ref($x) eq 'ARRAY';
  174   return $x if ref($x) eq $pkg;
  175   return $pkg->make(@{$x->data}) if Value::class($x) =~ m/Point|Vector/;
  176   Value::Error("Can't convert ".Value::showClass($x)." to a Matrix");
  177 }
  178 
  179 ############################################
  180 #
  181 #  Operations on matrices
  182 #
  183 
  184 sub add {
  185   my ($l,$r,$flag) = @_;
  186   if ($l->promotePrecedence($r)) {return $r->add($l,!$flag)}
  187   ($l,$r) = (promote($l)->data,promote($r)->data);
  188   Value::Error("Matrix addition with different dimensions")
  189     unless scalar(@{$l}) == scalar(@{$r});
  190   my @s = ();
  191   foreach my $i (0..scalar(@{$l})-1) {push(@s,$l->[$i] + $r->[$i])}
  192   return $pkg->make(@s);
  193 }
  194 
  195 sub sub {
  196   my ($l,$r,$flag) = @_;
  197   if ($l->promotePrecedence($r)) {return $r->sub($l,!$flag)}
  198   ($l,$r) = (promote($l)->data,promote($r)->data);
  199   Value::Error("Matrix subtraction with different dimensions")
  200     unless scalar(@{$l}) == scalar(@{$r});
  201   if ($flag) {my $tmp = $l; $l = $r; $r = $tmp};
  202   my @s = ();
  203   foreach my $i (0..scalar(@{$l})-1) {push(@s,$l->[$i] - $r->[$i])}
  204   return $pkg->make(@s);
  205 }
  206 
  207 sub mult {
  208   my ($l,$r,$flag) = @_;
  209   if ($l->promotePrecedence($r)) {return $r->mult($l,!$flag)}
  210   #
  211   #  Constant multiplication
  212   #
  213   if (Value::matchNumber($r) || Value::isComplex($r)) {
  214     my @coords = ();
  215     foreach my $x (@{$l->data}) {push(@coords,$x*$r)}
  216     return $pkg->make(@coords);
  217   }
  218   #
  219   #  Make points and vectors into columns if they are on the right
  220   #
  221   if (!$flag && Value::class($r) =~ m/Point|Vector/)
  222     {$r = (promote($r))->transpose} else {$r = promote($r)}
  223   #
  224   if ($flag) {my $tmp = $l; $l = $r; $r = $tmp}
  225   my @dl = $l->dimensions; my @dr = $r->dimensions;
  226   if (scalar(@dl) == 1) {@dl = (1,@dl); $l = $pkg->make($l)}
  227   if (scalar(@dr) == 1) {@dr = (@dr,1); $r = $pkg->make($r)->transpose}
  228   Value::Error("Can only multiply 2-dimensional matrices") if scalar(@dl) > 2 || scalar(@dr) > 2;
  229   Value::Error("Matices of dimensions $dl[0]x$dl[1] and $dr[0]x$dr[1] can't be multiplied")
  230     unless ($dl[1] == $dr[0]);
  231   #
  232   #  Do matrix multiplication
  233   #
  234   my @l = $l->value; my @r = $r->value;
  235   my @M = ();
  236   foreach my $j (0..$dr[1]-1) {
  237     my @row = ();
  238     foreach my $i (0..$dl[0]-1) {
  239       my $s = 0;
  240       foreach my $k (0..$dl[1]-1) {$s += $l[$i]->[$k] * $r[$k]->[$j]}
  241       push(@row,$s);
  242     }
  243     push(@M,$pkg->make(@row));
  244   }
  245   return $pkg->make(@M);
  246 }
  247 
  248 sub div {
  249   my ($l,$r,$flag) = @_;
  250   if ($l->promotePrecedence($r)) {return $r->div($l,!$flag)}
  251   Value::Error("Can't divide by a Matrix") if $flag;
  252   Value::Error("Matrices can only be divided by numbers")
  253     unless (Value::matchNumber($r) || Value::isComplex($r));
  254   Value::Error("Division by zero") if $r == 0;
  255   my @coords = ();
  256   foreach my $x (@{$l->data}) {push(@coords,$x/$r)}
  257   return $pkg->make(@coords);
  258 }
  259 
  260 sub power {
  261   my ($l,$r,$flag) = @_;
  262   if ($l->promotePrecedence($r)) {return $r->power($l,!$flag)}
  263   Value::Error("Can't use Matrices in exponents") if $flag;
  264   Value::Error("Only square matrices can be raised to a power") unless $l->isSquare;
  265   return Value::Matrix::I($l->length) if $r == 0;
  266   Value::Error("Matrix powers must be positive integers") unless $r =~ m/^[1-9]\d*$/;
  267   my $M = $l; foreach my $i (2..$r) {$M = $M*$l}
  268   return $M;
  269 }
  270 
  271 #
  272 #  Do lexicographic comparison
  273 #
  274 sub compare {
  275   my ($l,$r,$flag) = @_;
  276   if ($l->promotePrecedence($r)) {return $r->compare($l,!$flag)}
  277   ($l,$r) = (promote($l)->data,promote($r)->data);
  278   Value::Error("Matrix comparison with different dimensions")
  279     unless scalar(@{$l}) == scalar(@{$r});
  280   if ($flag) {my $tmp = $l; $l = $r; $r = $tmp};
  281   my $cmp = 0;
  282   foreach my $i (0..scalar(@{$l})-1) {
  283     $cmp = $l->[$i] <=> $r->[$i];
  284     last if $cmp;
  285   }
  286   return $cmp;
  287 }
  288 
  289 sub neg {
  290   my $p = promote(@_)->data;
  291   my @coords = ();
  292   foreach my $x (@{$p}) {push(@coords,-$x)}
  293   return $pkg->make(@coords);
  294 }
  295 
  296 #
  297 #  Transpose an  n x m  matrix
  298 #
  299 sub transpose {
  300   my $self = shift;
  301   my @d = $self->dimensions;
  302   if (scalar(@d) == 1) {@d = (1,@d); $self = $pkg->make($self)}
  303   Value::Error("Can't transpose ".scalar(@d)."-dimensional matrices") unless scalar(@d) == 2;
  304   my @M = (); my $M = $self->data;
  305   foreach my $j (0..$d[1]-1) {
  306     my @row = ();
  307     foreach my $i (0..$d[0]-1) {push(@row,$M->[$i]->data->[$j])}
  308     push(@M,$pkg->make(@row));
  309   }
  310   return $pkg->make(@M);
  311 }
  312 
  313 #
  314 #  Get an identity matrix of the requested size
  315 #
  316 sub I {
  317   my $d = shift; $d = shift if ref($d) eq $pkg;
  318   my @M = (); my @Z = split('',0 x $d);
  319   foreach my $i (0..$d-1) {
  320     my @row = @Z; $row[$i] = 1;
  321     push(@M,$pkg->make(@row));
  322   }
  323   return $pkg->make(@M);
  324 }
  325 
  326 #
  327 #  Extract a given row from the matrix
  328 #
  329 sub row {
  330   my $M = promote(shift); my $i = shift;
  331   return if $i == 0; $i-- if $i > 0;
  332   if ($M->isRow) {return if $i != 0; return $M}
  333   return $M->data->[$i];
  334 }
  335 
  336 #
  337 #  Extract a given element from the matrix
  338 #
  339 sub element {
  340   my $M = promote(shift);
  341   return $M->extract(@_);
  342 }
  343 
  344 #
  345 #  Extract a given column from the matrix
  346 #
  347 sub column {
  348   my $M = promote(shift); my $j = shift;
  349   return if $j == 0; $j-- if $j > 0;
  350   my @d = $M->dimensions; my @col = ();
  351   return if $j+1 > $d[1];
  352   return $M->data->[$j] if scalar(@d) == 1;
  353   foreach my $row (@{$M->data}) {push(@col,$pkg->make($row->data->[$j]))}
  354   return $pkg->make(@col);
  355 }
  356 
  357 # @@@ removeRow, removeColumn @@@
  358 # @@@ Det, inverse @@@
  359 
  360 ############################################
  361 #
  362 #  Generate the various output formats
  363 #
  364 
  365 sub stringify {
  366   my $self = shift;
  367   return $self->TeX(undef,$self->{open},$self->{close}) if $$Value::context->flag('StringifyAsTeX');
  368   return $self->string(undef,$self->{open},$self->{close})
  369 }
  370 
  371 sub string {
  372   my $self = shift; my $equation = shift;
  373   my $def = ($equation->{context} || $$Value::context)->lists->get('Matrix');
  374   my $open = shift || $def->{open}; my $close = shift || $def->{close};
  375   my @coords = ();
  376   foreach my $x (@{$self->data}) {
  377     if (Value::isValue($x)) {push(@coords,$x->string($equation,$open,$close))}
  378       else {push(@coords,$x)}
  379   }
  380   return $open.join(',',@coords).$close;
  381 }
  382 
  383 #
  384 #  Use \matrix to lay out matrices
  385 #
  386 sub TeX {
  387   my $self = shift; my $equation = shift;
  388   my $def = ($equation->{context} || $$Value::context)->lists->get('Matrix');
  389   my $open = shift || $def->{open}; my $close = shift || $def->{close};
  390   $open = '\{' if $open eq '{'; $close = '\}' if $close eq '}';
  391   my $TeX = ''; my @entries = (); my $d;
  392   if ($self->isRow) {
  393     foreach my $x (@{$self->data}) {
  394       push(@entries,(Value::isValue($x))? $x->TeX($equation): $x);
  395     }
  396     $TeX .= join(' &',@entries) . "\n";
  397     $d = scalar(@entries);
  398   } else {
  399     foreach my $row (@{$self->data}) {
  400       foreach my $x (@{$row->data}) {
  401         push(@entries,(Value::isValue($x))? $x->TeX($equation): $x);
  402       }
  403       $TeX .= join(' &',@entries) . '\cr'."\n";
  404       $d = scalar(@entries); @entries = ();
  405     }
  406   }
  407   return '\left'.$open.'\begin{array}{'.('c'x$d).'}'."\n".$TeX.'\end{array}\right'.$close;
  408 }
  409 
  410 ###########################################################################
  411 
  412 1;
  413 

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9