[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 6899 - (download) (as text) (annotate)
Sat Jun 25 20:51:08 2011 UTC (8 years, 5 months ago) by dpvc
File size: 15940 byte(s)
Initial attempt to link Matrix MathObject to CPAN Matrix library (added methods for the major ones in the CPAN library)

    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; no strict "refs";
   11 use Matrix; use Complex1;
   12 our @ISA = qw(Value);
   13 
   14 #
   15 #  Convert a value to a matrix.  The value can be:
   16 #     a list of numbers or list of (nested) references to arrays of numbers,
   17 #     a point, vector or matrix object, a matrix-valued formula, or a string
   18 #     that evaluates to a matrix
   19 #
   20 sub new {
   21   my $self = shift; my $class = ref($self) || $self;
   22   my $context = (Value::isContext($_[0]) ? shift : $self->context);
   23   my $M = shift; $M = [] unless defined $M; $M = [$M,@_] if scalar(@_) > 0;
   24   $M = @{$M}[0] if ref($M) =~ m/^Matrix(Real1)?/;
   25   $M = Value::makeValue($M,context=>$context) if ref($M) ne 'ARRAY';
   26   return bless {data => $M->data, context=>$context}, $class
   27     if (Value::classMatch($M,'Point','Vector','Matrix') && scalar(@_) == 0);
   28   return $M if Value::isFormula($M) && Value::classMatch($self,$M->type);
   29   my @M = (ref($M) eq 'ARRAY' ? @{$M} : $M);
   30   Value::Error("Matrices must have at least one entry") unless scalar(@M) > 0;
   31   return $self->matrixMatrix($context,@M) if ref($M[0]) eq 'ARRAY' ||
   32     Value::classMatch($M[0],'Matrix','Vector','Point') ||
   33     (Value::isFormula($M[0]) && $M[0]->type =~ m/Matrix|Vector|Point/);
   34   return $self->numberMatrix($context,@M);
   35 }
   36 
   37 #
   38 #  (Recursively) make a matrix from a list of array refs
   39 #  and report errors about the entry types
   40 #
   41 sub matrixMatrix {
   42   my $self = shift; my $class = ref($self) || $self;
   43   my $context = shift;
   44   my ($x,$m); my @M = (); my $isFormula = 0;
   45   foreach $x (@_) {
   46     if (Value::isFormula($x)) {push(@M,$x); $isFormula = 1} else {
   47       $m = $self->new($context,$x); push(@M,$m);
   48       $isFormula = 1 if Value::isFormula($m);
   49     }
   50   }
   51   my ($type,$len) = ($M[0]->entryType->{name},$M[0]->length);
   52   foreach $x (@M) {
   53     Value::Error("Matrix rows must all be the same type")
   54       unless (defined($x->entryType) && $type eq $x->entryType->{name});
   55     Value::Error("Matrix rows must all be the same length") unless ($len eq $x->length);
   56   }
   57   return $self->formula([@M]) if $isFormula;
   58   bless {data => [@M], context=>$context}, $class;
   59 }
   60 
   61 #
   62 #  Form a 1 x n matrix from a list of numbers
   63 #  (could become a row of an  m x n  matrix)
   64 #
   65 sub numberMatrix {
   66   my $self = shift; my $class = ref($self) || $self;
   67   my $context = shift;
   68   my @M = (); my $isFormula = 0;
   69   foreach my $x (@_) {
   70     $x = Value::makeValue($x,context=>$context);
   71     Value::Error("Matrix row entries must be numbers") unless Value::isNumber($x);
   72     push(@M,$x); $isFormula = 1 if Value::isFormula($x);
   73   }
   74   return $self->formula([@M]) if $isFormula;
   75   bless {data => [@M], context=>$context}, $class;
   76 }
   77 
   78 #
   79 #  Recursively get the entries in the matrix and return
   80 #  an array of (references to arrays of ... ) numbers
   81 #
   82 sub value {
   83   my $self = shift;
   84   my $M = $self->data;
   85   return @{$M} unless Value::classMatch($M->[0],'Matrix');
   86   my @M = ();
   87   foreach my $x (@{$M}) {push(@M,[$x->value])}
   88   return @M;
   89 }
   90 
   91 #
   92 #  Recursively get the dimensions of the matrix.
   93 #  Returns (n) for a 1 x n, or (n,m) for an n x m, etc.
   94 #
   95 sub dimensions {
   96   my $self = shift;
   97   my $r = $self->length;
   98   my $v = $self->data;
   99   return ($r,) unless Value::classMatch($v->[0],'Matrix');
  100   return ($r,$v->[0]->dimensions);
  101 }
  102 
  103 #
  104 #  Return the proper type for the matrix
  105 #
  106 sub typeRef {
  107   my $self = shift;
  108   return Value::Type($self->class, $self->length, $Value::Type{number})
  109     unless Value::classMatch($self->data->[0],'Matrix');
  110   return Value::Type($self->class, $self->length, $self->data->[0]->typeRef);
  111 }
  112 
  113 #
  114 #  True if the matrix is a square matrix
  115 #
  116 sub isSquare {
  117   my $self = shift;
  118   my @d = $self->dimensions;
  119   return 0 if scalar(@d) > 2;
  120   return 1 if scalar(@d) == 1 && $d[0] == 1;
  121   return $d[0] == $d[1];
  122 }
  123 
  124 #
  125 #  True if the matrix is 1-dimensional (i.e., is a matrix row)
  126 #
  127 sub isRow {
  128   my $self = shift;
  129   my @d = $self->dimensions;
  130   return scalar(@d) == 1;
  131 }
  132 
  133 #
  134 #  See if the matrix is an Indenity matrix
  135 #
  136 sub isOne {
  137   my $self = shift;
  138   return 0 unless $self->isSquare;
  139   my $i = 0;
  140   foreach my $row (@{$self->{data}}) {
  141     my $j = 0;
  142     foreach my $k (@{$row->{data}}) {
  143       return 0 unless $k eq (($i == $j)? "1": "0");
  144       $j++;
  145     }
  146     $i++;
  147   }
  148   return 1;
  149 }
  150 
  151 #
  152 #  See if the matrix is all zeros
  153 #
  154 sub isZero {
  155   my $self = shift;
  156   foreach my $x (@{$self->{data}}) {return 0 unless $x->isZero}
  157   return 1;
  158 }
  159 
  160 #
  161 #  Make arbitrary data into a matrix, if possible
  162 #
  163 sub promote {
  164   my $self = shift; my $class = ref($self) || $self;
  165   my $context = (Value::isContext($_[0]) ? shift : $self->context);
  166   my $x = (scalar(@_) ? shift : $self);
  167   return $self->new($context,$x,@_) if scalar(@_) > 0 || ref($x) eq 'ARRAY';
  168   $x = Value::makeValue($x,context=>$context);
  169   return $x->inContext($context) if ref($x) eq $class;
  170   return $self->make($context,@{$x->data}) if Value::classMatch($x,'Point','Vector');
  171   Value::Error("Can't convert %s to %s",Value::showClass($x),Value::showClass($self));
  172 }
  173 
  174 #
  175 #  Don't inherit ColumnVector flag
  176 #
  177 sub noinherit {
  178   my $self = shift;
  179   return ("ColumnVector",$self->SUPER::noinherit);
  180 }
  181 
  182 ############################################
  183 #
  184 #  Operations on matrices
  185 #
  186 
  187 sub add {
  188   my ($self,$l,$r,$other) = Value::checkOpOrderWithPromote(@_);
  189   my @l = @{$l->data}; my @r = @{$r->data};
  190   Value::Error("Can't add Matrices with different dimensions")
  191     unless scalar(@l) == scalar(@r);
  192   my @s = ();
  193   foreach my $i (0..scalar(@l)-1) {push(@s,$l[$i] + $r[$i])}
  194   return $self->inherit($other)->make(@s);
  195 }
  196 
  197 sub sub {
  198   my ($self,$l,$r,$other) = Value::checkOpOrderWithPromote(@_);
  199   my @l = @{$l->data}; my @r = @{$r->data};
  200   Value::Error("Can't subtract Matrices with different dimensions")
  201     unless scalar(@l) == scalar(@r);
  202   my @s = ();
  203   foreach my $i (0..scalar(@l)-1) {push(@s,$l[$i] - $r[$i])}
  204   return $self->inherit($other)->make(@s);
  205 }
  206 
  207 sub mult {
  208   my ($l,$r,$flag) = @_; my $self = $l; my $other = $r;
  209   #
  210   #  Constant multiplication
  211   #
  212   if (Value::matchNumber($r) || Value::isComplex($r)) {
  213     my @coords = ();
  214     foreach my $x (@{$l->data}) {push(@coords,$x*$r)}
  215     return $self->make(@coords);
  216   }
  217   #
  218   #  Make points and vectors into columns if they are on the right
  219   #
  220   if (!$flag && Value::classMatch($r,'Point','Vector'))
  221     {$r = ($self->promote($r))->transpose} else {$r = $self->promote($r)}
  222   #
  223   if ($flag) {my $tmp = $l; $l = $r; $r = $tmp}
  224   my @dl = $l->dimensions; my @dr = $r->dimensions;
  225   if (scalar(@dl) == 1) {@dl = (1,@dl); $l = $self->make($l)}
  226   if (scalar(@dr) == 1) {@dr = (@dr,1); $r = $self->make($r)->transpose}
  227   Value::Error("Can only multiply 2-dimensional matrices") if scalar(@dl) > 2 || scalar(@dr) > 2;
  228   Value::Error("Matices of dimensions %dx%d and %dx%d can't be multiplied",@dl,@dr)
  229     unless ($dl[1] == $dr[0]);
  230   #
  231   #  Do matrix multiplication
  232   #
  233   my @l = $l->value; my @r = $r->value; my @M = ();
  234   foreach my $i (0..$dl[0]-1) {
  235     my @row = ();
  236     foreach my $j (0..$dr[1]-1) {
  237       my $s = 0;
  238       foreach my $k (0..$dl[1]-1) {$s += $l[$i]->[$k] * $r[$k]->[$j]}
  239       push(@row,$s);
  240     }
  241     push(@M,$self->make(@row));
  242   }
  243   $self = $self->inherit($other) if Value::isValue($other);
  244   return $self->make(@M);
  245 }
  246 
  247 sub div {
  248   my ($l,$r,$flag) = @_; my $self = $l;
  249   Value::Error("Can't divide by a Matrix") if $flag;
  250   Value::Error("Matrices can only be divided by Numbers")
  251     unless (Value::matchNumber($r) || Value::isComplex($r));
  252   Value::Error("Division by zero") if $r == 0;
  253   my @coords = ();
  254   foreach my $x (@{$l->data}) {push(@coords,$x/$r)}
  255   return $self->make(@coords);
  256 }
  257 
  258 sub power {
  259   my ($l,$r,$flag) = @_; my $self = shift; my $context = $self->context;
  260   Value::Error("Can't use Matrices in exponents") if $flag;
  261   Value::Error("Only square matrices can be raised to a power") unless $l->isSquare;
  262   $r = Value::makeValue($r,context=>$context);
  263   Value::Error("Matrix powers must be non-negative integers") unless $r->isNumber && $r =~ m/^\d+$/;
  264   return $context->Package("Matrix")->I($l->length,$context) if $r == 0;
  265   my $M = $l; foreach my $i (2..$r) {$M = $M*$l}
  266   return $M;
  267 }
  268 
  269 #
  270 #  Do lexicographic comparison (row by row)
  271 #
  272 sub compare {
  273   my ($self,$l,$r) = Value::checkOpOrderWithPromote(@_);
  274   Value::Error("Can't compare Matrices with different dimensions")
  275     unless join(',',$l->dimensions) eq join(',',$r->dimensions);
  276   my @l = @{$l->data}; my @r = @{$r->data};
  277   foreach my $i (0..scalar(@l)-1) {
  278     my $cmp = $l[$i] <=> $r[$i];
  279     return $cmp if $cmp;
  280   }
  281   return 0;
  282 }
  283 
  284 sub neg {
  285   my $self = promote(@_); my @coords = ();
  286   foreach my $x (@{$self->data}) {push(@coords,-$x)}
  287   return $self->make(@coords);
  288 }
  289 
  290 #
  291 #  Transpose an  n x m  matrix
  292 #
  293 sub transpose {
  294   my $self = promote(@_);
  295   my @d = $self->dimensions;
  296   if (scalar(@d) == 1) {@d = (1,@d); $self = $self->make($self)}
  297   Value::Error("Can't transpose %d-dimensional matrices",scalar(@d)) unless scalar(@d) == 2;
  298   my @M = (); my $M = $self->data;
  299   foreach my $j (0..$d[1]-1) {
  300     my @row = ();
  301     foreach my $i (0..$d[0]-1) {push(@row,$M->[$i]->data->[$j])}
  302     push(@M,$self->make(@row));
  303   }
  304   return $self->make(@M);
  305 }
  306 
  307 #
  308 #  Get an identity matrix of the requested size
  309 #
  310 sub I {
  311   my $self = shift; my $d = shift; my $context = shift || $self->context;
  312   $d = ($self->dimensions)[0] if !defined $d && ref($self) && $self->isSquare;
  313   Value::Error("You must provide a dimension for the Identity matrix") unless defined $d;
  314   Value::Error("Dimension must be a positive integer") unless $d =~ m/^[1-9]\d*$/;
  315   my @M = (); my @Z = split('',0 x $d);
  316   foreach my $i (0..$d-1) {
  317     my @row = @Z; $row[$i] = 1;
  318     push(@M,$self->make($context,@row));
  319   }
  320   return $self->make($context,@M);
  321 }
  322 
  323 #
  324 #  Extract a given row from the matrix
  325 #
  326 sub row {
  327   my $self = (ref($_[0]) ? $_[0] : shift);
  328   my $M = $self->promote(shift); my $i = shift;
  329   return if $i == 0; $i-- if $i > 0;
  330   if ($M->isRow) {return if $i != 0 && $i != -1; return $M}
  331   return $M->data->[$i];
  332 }
  333 
  334 #
  335 #  Extract a given column from the matrix
  336 #
  337 sub column {
  338   my $self = (ref($_[0]) ? $_[0] : shift);
  339   my $M = $self->promote(shift); my $j = shift;
  340   return if $j == 0; $j-- if $j > 0;
  341   my @d = $M->dimensions;
  342   if (scalar(@d) == 1) {
  343     return if $j+1 > $d[0] || $j < -$d[0];
  344     return $M->data->[$j];
  345   }
  346   return if $j+1 > $d[1] || $j < -$d[1];
  347   my @col = ();
  348   foreach my $row (@{$M->data}) {push(@col,$self->make($row->data->[$j]))}
  349   return $self->make(@col);
  350 }
  351 
  352 #
  353 #  Extract a given element from the matrix
  354 #
  355 sub element {
  356   my $self = (ref($_[0]) ? $_[0] : shift);
  357   my $M = $self->promote(shift);
  358   return $M->extract(@_);
  359 }
  360 
  361 # @@@ assign @@@
  362 # @@@ removeRow, removeColumn @@@
  363 # @@@ Minor @@@
  364 
  365 
  366 ##################################################
  367 #
  368 #  Convert MathObject Matrix to old-style Matrix
  369 #
  370 sub wwMatrix {
  371   my $self = (ref($_[0]) ? $_[0] : shift);
  372   my $M = $self->promote(shift); my $j = shift; my $wwM;
  373   return $self->{wwM} if defined($self->{wwM});
  374   my @d = $M->dimensions;
  375   Value->Error("Matrix must be two-dimensional to convert to MatrixReal1") if scalar(@d) > 2;
  376   if (scalar(@d) == 1) {
  377     $wwM = new Matrix(1,$d[0]);
  378     foreach my $j (0..$d[0]-1) {
  379       $wwM->[0][0][$j] = $self->wwMatrixEntry($M->data->[$j]);
  380     }
  381   } else {
  382     $wwM = new Matrix(@d);
  383     foreach my $i (0..$d[0]-1) {
  384       my $row = $M->data->[$i];
  385       foreach my $j (0..$d[1]-1) {
  386   $wwM->[0][$i][$j] = $self->wwMatrixEntry($row->data->[$j]);
  387       }
  388     }
  389   }
  390   $self->{wwM} = $wwM;
  391   return $wwM;
  392 }
  393 sub wwMatrixEntry {
  394   my $self = shift; my $x = shift;
  395   return $x->value if $x->isReal;
  396   return Complex1::cplx($x->Re->value,$x->Im->value) if $x->isComplex;
  397   return $x;
  398 }
  399 sub wwMatrixLR {
  400   my $self = shift;
  401   return $self->{lrM} if defined($self->{lrM});
  402   $self->wwMatrix;
  403   $self->{lrM} = $self->{wwM}->decompose_LR;
  404   return $self->{lrM};
  405 }
  406 
  407 
  408 ###################################
  409 #
  410 #  From MatrixReal1.pm
  411 #
  412 
  413 sub det {
  414   my $self = shift; $self->wwMatrixLR;
  415   return $self->{lrM}->det_LR;
  416 }
  417 
  418 sub inverse {
  419   my $self = shift; $self->wwMatrixLR;
  420   return $self->new($self->{lrM}->invert_LR);
  421 }
  422 
  423 sub decompose_LR {
  424   my $self = shift;
  425   return $self->wwMatrixLR;
  426 }
  427 
  428 sub dim {
  429   my $self = shift;
  430   return $self->wwMatrix->dim();
  431 }
  432 
  433 sub norm_one {
  434   my $self = shift;
  435   return $self->wwMatrix->norm_one();
  436 }
  437 
  438 sub norm_max {
  439   my $self = shift;
  440   return $self->wwMatrix->norm_max();
  441 }
  442 
  443 sub kleene {
  444   my $self = shift;
  445   return $self->new($self->wwMatrix->kleene());
  446 }
  447 
  448 sub normalize {
  449   my $self = shift;
  450   my $v = $self->new(shift)->wwMatrix;
  451   my ($M,$b) = $self->wwMatrix->normalize($v);
  452   return ($self->new($M),$self->new($b));
  453 }
  454 
  455 sub solve_LR {
  456   my $self = shift;
  457   my $v = $self->new(shift)->wwMatrix;
  458   my ($d,$b,$M) = $self->lrMatrix->solve_LR($v);
  459   return ($d,$self->new($b),$self->new($M));
  460 }
  461 
  462 sub condition {
  463   my $self = shift;
  464   my $I = $self->new(shift)->wwMatrix;
  465   return $self->new($self->wwMatrix->condition($I));
  466 }
  467 
  468 sub order_LR {
  469   my $self = shift;
  470   return $self->wwMatrixLR->order_LR;
  471 }
  472 
  473 sub solve_GSM {
  474   my $self = shift;
  475   my $x0 = $self->new(shift)->wwMatrix;
  476   my $b  = $self->new(shift)->wwMatrix;
  477   my $e  = shift;
  478   my $v = $self->wwMatrix->solve_GSM($x0,$b,$e);
  479   $v = $self->new($v) if defined($v);
  480   return $v;
  481 }
  482 
  483 sub solve_SSM {
  484   my $self = shift;
  485   my $x0 = $self->new(shift)->wwMatrix;
  486   my $b  = $self->new(shift)->wwMatrix;
  487   my $e  = shift;
  488   my $v = $self->wwMatrix->solve_SSM($x0,$b,$e);
  489   $v = $self->new($v) if defined($v);
  490   return $v;
  491 }
  492 
  493 sub solve_RM {
  494   my $self = shift;
  495   my $x0 = $self->new(shift)->wwMatrix;
  496   my $b  = $self->new(shift)->wwMatrix;
  497   my $w = shift; my $e = shift;
  498   my $v = $self->wwMatrix->solve_RM($x0,$b,$w,$e);
  499   $v = $self->new($v) if defined($v);
  500   return $v;
  501 }
  502 
  503 sub is_symmetric {
  504   my $self = shift;
  505   return $self->wwMatrix->is_symmetric;
  506 }
  507 
  508 ###################################
  509 #
  510 #  From Matrix.pm
  511 #
  512 
  513 sub trace {
  514   my $self = shift;
  515   return $self->wwMatrix->trace;
  516 }
  517 
  518 sub proj {
  519   my $self = shift;
  520   my $v = $self->new(shift)->wwMatrix;
  521   return $self->new($self->wwMatrix->proj($v));
  522 }
  523 
  524 sub proj_coeff {
  525   my $self = shift;
  526   my $v = $self->new(shift)->wwMatrix;
  527   return $self->new($self->wwMatrix->proj_coeff($v));
  528 }
  529 
  530 sub L {
  531   my $self = shift;
  532   return $self->new($self->wwMatrixLR->L);
  533 }
  534 
  535 sub R {
  536   my $self = shift;
  537   return $self->new($self->wwMatrixLR->R);
  538 }
  539 
  540 sub PL {
  541   my $self = shift;
  542   return $self->new($self->wwMatrixLR->PL);
  543 }
  544 
  545 sub PR {
  546   my $self = shift;
  547   return $self->new($self->wwMatrixLR->PR);
  548 }
  549 
  550 
  551 ############################################
  552 #
  553 #  Generate the various output formats
  554 #
  555 
  556 #
  557 #  Use array environment to lay out matrices
  558 #
  559 sub TeX {
  560   my $self = shift; my $equation = shift;
  561   my $def = ($equation->{context} || $self->context)->lists->get('Matrix');
  562   my $open  = shift || $self->{open} || $def->{open};
  563   my $close = shift || $self->{close} || $def->{close};
  564   $open =~ s/([{}])/\\$1/g; $close =~ s/([{}])/\\$1/g;
  565   my $TeX = ''; my @entries = (); my $d;
  566   if ($self->isRow) {
  567     foreach my $x (@{$self->data}) {
  568       if (Value::isValue($x)) {
  569   $x->{format} = $self->{format} if defined $self->{format};
  570   push(@entries,$x->TeX($equation));
  571       } else {push(@entries,$x)}
  572     }
  573     $TeX .= join(' &',@entries) . "\n";
  574     $d = scalar(@entries);
  575   } else {
  576     foreach my $row (@{$self->data}) {
  577       foreach my $x (@{$row->data}) {
  578   if (Value::isValue($x)) {
  579     $x->{format} = $self->{format} if defined $self->{format};
  580     push(@entries,$x->TeX($equation));
  581   } else {push(@entries,$x)}
  582       }
  583       $TeX .= join(' &',@entries) . '\cr'."\n";
  584       $d = scalar(@entries); @entries = ();
  585     }
  586   }
  587   $TeX =~ s/\\cr\n$/\n/;
  588   return '\left'.$open.'\begin{array}{'.('c'x$d).'}'."\n".$TeX.'\end{array}\right'.$close;
  589 }
  590 
  591 ###########################################################################
  592 
  593 1;
  594 

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9