3
PDL::Matrix -- a derived matrix class that implements column-major
4
constructors and methods
8
This document refers to version PDL::Matrix 0.01 of PDL::Matrix
14
$m = mpdl [[1,2,3],[4,5,6]];
15
$m = PDL::Matrix->pdl([[1,2,3],[4,5,6]]);
17
@dimsa = $a->mdims; # 'dims' is not overloaded
26
This package tries to help people who want to use PDL for 2D matrix
27
computation with lots of indexing involved . It provides a PDL
28
subclass so one- and two-dimensional piddles that are used as
29
vectors resp. matrices can be typed in using traditional matrix
32
The original pdl class refers to the first index as the first row,
33
the second index as the first column of a matrix. Consider
35
print $B = sequence(3,2)
41
which gives a 2x3 matrix in terms of the matrix convention, but the
42
constructor used (3,2). This might get more confusing when using
43
slices like sequence(3,2)->slice("1:2,(0)") : with traditional
44
matrix convention one would expect [2 4] instead of [1 2].
46
This subclass PDL::Matrix overloads the constructors and indexing
47
functions of pdls so that they are compatible with the usual matrix
48
convention, where the first dimension refers to the row of a
49
matrix. So now, the above example would be written as
51
print $B = PDL::Matrix->sequence(3,2) # or $B = msequence(3,2)
58
Routines like eigenvalue or matrix inversion can be used without
61
Furthermore one can construct and use vectors as n x 1 matrices
62
without mentioning the second index '1'.
66
C<PDL::Matrix> works by overloading a number of PDL constructors
67
and methods such that first and second args (corresponding to
68
first and second dims of corresponding matrices) are effectively swapped.
69
It is not yet clear if PDL::Matrix achieves a consistent column
70
major look-and-feel in this way.
85
@ISA = qw/PDL::Exporter PDL/;
89
#######################################################################=
94
# --------> constructors
96
=head2 mpdl, PDL::Matrix::pdl
100
constructs an object of class PDL::Matrix which is a piddle child
101
class, where the first index refers to the first column of the
102
two-dimensional piddle.
106
$m = mpdl [[1,2,3],[4,5,6]];
107
$m = PDL::Matrix->pdl([[1,2,3],[4,5,6]]);
113
my $pdl = $class->SUPER::pdl(@_);
114
bless $pdl, ref $class || $class;
117
=head2 mzeroes, mones, msequence
121
constructs a PDL::Matrix object similar to the piddle constructors
122
zeroes, ones, sequence
126
# for constructors with specified dimensions I only have to swap these
128
for my $func (qw /zeroes ones sequence/) {
134
ref(\$arg[0]) ne 'PDL::Type' ? (\@arg > 1 ? (\@arg[1,0] = \@arg[0,1]) :
135
(\@arg[0,1] = (1,\$arg[0])) ) :
136
(\@arg > 2 ? (\@arg[2,1] = \@arg[1,2]) :
137
(\@arg[1,2] = (1,\$arg[1])) );
138
my \$pdl = \$class->SUPER::$func(\@arg);
139
bless \$pdl, ref \$class || \$class;
143
# print "evaluating $code\n";
147
# functions that construct a matrix pdl and that can be exported, very
148
# trivial, they just call its methods optional type argument is
149
# checked within these methods
150
for my $func (qw /pdl zeroes ones sequence/) {
151
push @EXPORT_OK, "m$func";
154
sub m$func { PDL::Matrix->$func(\@_) }
157
# print "evaluating $code\n";
165
constructs an object of class PDL::Matrix which is of matrix
170
print $v = vpdl [0,1];
179
my $pdl = transpose(PDL->pdl(@_));
180
bless $pdl, PDL::Matrix;
182
push @EXPORT_OK, "vpdl";
184
=head2 vzeroes, vones, vsequence
188
constructs a PDL::Matrix object with matrix dimensions (n x 1),
189
therefore only the first scalar argument is used.
193
print $v = vsequence(2);
201
for my $func (qw /zeroes ones sequence/) {
202
push @EXPORT_OK, "v$func";
207
ref(\$arg[0]) ne 'PDL::Type' ? (\@arg = (\$arg[0],1)) :
208
(\@arg = (\$arg[0],\$arg[1],1));
209
PDL::Matrix->$func(\@arg);
213
# print "evaluating $code\n";
219
=head2 PDL::Matrix::slice, PDL::Matrix::dice
223
same as slice, dice for normal piddles, but reflecting the matrix
224
convention by swapping the first two arguments.
228
print sequence(3,2)->slice("1:2,(0)") # piddle
230
print msequence(3,2)->slice("1:2,(0)") # PDL::Matrix
238
# add another dimension if slices has only one
239
# convenient for vectors
240
$ind =~ s/^([^,]*)$/$1,/;
241
# swap first two arguments
242
$ind =~ s/^([^,]*),([^,]*)(.*)/$2,$1$3/;
243
$self->SUPER::slice($ind);
251
@arg >=2 ? @arg[1,0] = @arg[0,1] : ( (@arg == 1 && $self->dims == 2) ?
252
@arg = (0,$arg[0]) : 1 );
253
$self->SUPER::dice(@arg);
257
=head2 PDL::Matrix::at
261
same as at for piddles, but reflecting the matrix convention by
262
swapping the first two arguments
264
If only one scalar argument is used, we assume the object to be a
265
vector and look only at the first column.
269
# swap arguments if number of arguments is greater than 1
270
# if its one, look at the first column, assuming it is a vector
274
@arg >=2 ? @arg[1,0] = @arg[0,1] : ( (@arg == 1 && $self->dims == 2) ?
275
@arg = (0,$arg[0]) : 1 );
276
$self->SUPER::at(@arg);
279
=head2 PDL::Matrix::set
283
set a particular value in a PDL::Matrix object. Note that this has to
284
be called as an object method rather than a function
288
print msequence(3,3)->set(2,0,-1) # ok with PDL::Matrix convention
295
print set msequence(3,3), 2,0,-1 # does not conform with PDL::Matrix convention
305
# this does not work as the usual set, but can be used via
306
# $M->set(0,2,0.4) which sets the element in the first row and third
308
sub set{ # Sets a particular single value
309
barf("Usage: set($pdl, $x, $y,.., $value)") if $#_<2;
310
my $self = shift; my $value = pop @_;
312
@arg >=2 ? @arg[1,0] = @arg[0,1] : ( (@arg == 1 && $self->dims == 2) ?
313
@arg = (0,$arg[0]) : 1 );
314
PDL::Core::set_c ($self, [@arg], $value);
318
=head2 PDL::Matrix::reshape
322
same as reshape for piddles, but reflecting the matrix convention by
323
swapping the first two arguments
330
@arg >=2 ? @arg[1,0] = @arg[0,1] : (@arg = (1,$arg[0]));
331
$self->SUPER::reshape(@arg);
334
# this is needed because only doing
335
# > use overload '~' => \&PDL::transpose;
336
# would not give an object of type PDL::Matrix back
337
# => possible bug in transpose?!
340
my $pdl = $self->PDL::transpose;
341
bless $pdl, ref $self;
344
use overload '~' => \&PDL::Matrix::transpose;
350
returns the dimensions of the PDL::Matrix object in matrix
353
C<dims> is NOT overloaded by PDL::Matrix to make sure that
354
methods like PDL::transpose still work. So use C<mdims> to get
355
the dims in the PDL::Matrix notation.
359
print msequence(3,2)->mdims
364
# I cannot overload dims, because it is apparently used many times
365
# in methods like transpose! :-(
366
# this is not nice but I don't know how to avoid this
369
my @res = $self->SUPER::dims;
370
@res >=2 ? @res[1,0] = @res[0,1] : 1;
376
return matinv($self);
380
=head2 kroneckerproduct
384
returns kroneckerproduct of two matrices. This is not efficiently
388
print kroneckerproduct(msequence(2,2),mones(2,2))
398
# returns kroneckerproduct of two matrices
399
sub kroneckerproduct {
402
my ($r0,$c0) = $arg[0]->mdims;
403
my ($r1,$c1) = $arg[1]->mdims;
405
my $out = mzeroes($r0*$r1,$c0*$c1);
407
for (my $i=0;$i<$r0;$i++) {
408
for (my $j=0;$j<$c0;$j++) {
409
($_ = $out->slice(($i*$r1).":".(($i+1)*$r1-1).",".
410
($j*$c1).":".(($j+1)*$c1-1)) ) .= $arg[0]->at($i,$j) * $arg[1];
416
push @EXPORT_OK, "kroneckerproduct";
419
my ($self,@args) = @_;
420
return $self->transpose->SUPER::rotate(@args)->transpose;
426
return PDL::sumover(transpose($mpdl)->xchg(0,2));
428
push @EXPORT_OK, "msumover";
435
returns a generalized determinant of a matrix. If the matrix is not
436
regular, one can specify the rank of the matrix and the corresponding
437
subdeterminant is returned. This is implemented using the C<eigens>
441
print msequence(3,3)->determinant(2) # determinant of
442
# regular 2x2 submatrix
449
my ($mpdl,$rank) = @_;
450
my $eigenvalues = (PDL::Math::eigens($mpdl))[1];
451
my @sort = list(PDL::Ufunc::qsorti(abs($eigenvalues)));
452
$eigenvalues = $eigenvalues->dice([@sort[-$rank..-1]]);
453
PDL::Ufunc::dprod($eigenvalues);
460
returns the trace of a matrix (sum of diagonals)
466
$mpdl->diagonal(0,1)->sum;
469
# this has to be overloaded so that the PDL::slice
470
# is called and not PDL::Matrix::slice :-(
473
$dim = $pdl->getndims+1+$dim if $dim < 0;
474
barf ("too high/low dimension in call to dummy, allowed min/max=0/"
476
if $dim>$pdl->getndims || $dim < 0;
477
$_[2] = 1 if ($#_ < 2);
478
$pdl->PDL::slice((','x$dim)."*$_[2]");
482
# now some of my very own helper functions...
483
# stupid function to print a PDL::Matrix object in Maple code
485
my ($self,@args) = @_;
487
my ($dimR,$dimC) = mdims($self);
490
$s .= $args[0].":=" unless $args[0] eq "";
491
if (defined($dimR)) {
492
$s .= "matrix($dimR,$dimC,[";
493
for(my $i=0;$i<$dimR;++$i) {
495
for(my $j=0;$j<$dimC;++$j) {
496
$s .= $self->at($i,$j);
497
$s .= "," if $j+1<$dimC;
500
$s .= "," if $i+1<$dimR;
505
$s = "vector($dimC,[";
506
for(my $i=0;$i<$dimC;++$i) {
508
$s .= "," if $i+1<$dimC;
515
print stringifymaple(@_).";\n";
518
# stupid function to print a PDL::Matrix object in (La)TeX code
520
my ($self,@args) = @_;
522
my ($dimR,$dimC) = mdims($self);
525
$s .= $args[0]."=" unless $args[0] eq "";
526
$s .= "\\begin{pmatrix}\n";
527
for(my $i=0;$i<$dimR;++$i) {
528
for(my $j=0;$j<$dimC;++$j) {
529
$s .= $self->at($i,$j);
530
$s .= " & " if $j+1<$dimC;
532
$s .= " \\\\ \n" if $i+1<$dimR;
534
$s .= "\n \\end{pmatrix}\n";
540
print stringifyTeX(@_)."\n";
544
=head2 vcrossp, PDL::Matrix::crossp
548
similar to PDL::crossp, however reflecting PDL::Matrix notations
552
# crossp for my special vectors
554
my ($pdl1,$pdl2) = @_;
555
return transpose(PDL::crossp(~$pdl1,~$pdl2));
557
sub vcrossp { PDL::Matrix->crossp(\@_) }
558
push @EXPORT_OK, "vcrossp";
560
%EXPORT_TAGS = (Func=>[@EXPORT_OK]);
564
=head1 BUGS AND PROBLEMS
566
Because we change the way piddles are constructed, not all pdl
567
operators may be applied to piddle-matrices. The inner product is not
568
redefined. We might have missed some functions/methods. Internal
569
consistency of our approach needs yet to be established.
573
check all PDL functions, benchmarks, optimization, lots of other things ...
577
Stephan Heuel (stephan@heuel.org), Christian Soeller
578
(c.soeller@auckland.ac.nz).
582
All rights reserved. There is no warranty. You are allowed to
583
redistribute this software / documentation under certain
584
conditions. For details, see the file COPYING in the PDL
585
distribution. If this file is separated from the PDL distribution, the
586
copyright notice should be included in the file.
added
removed
Basic/Matrix.pm
