[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 5042 - (download) (as text) (annotate)
Thu Jun 28 01:31:09 2007 UTC (12 years, 8 months ago) by dpvc
File size: 11567 byte(s)
Recent changes to automatically do promotion in the Value methods was
a mistake.  I put it back into the subclass methods again.

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9