File Coverage

blib/lib/PDL/Matrix.pm
Criterion Covered Total %
statement 19 102 18.6
branch 3 38 7.8
condition 1 6 16.6
subroutine 7 27 25.9
pod 12 22 54.5
total 42 195 21.5


line stmt bran cond sub pod time code
1             =head1 NAME
2              
3             PDL::Matrix -- a convenience matrix class for column-major access
4              
5             =head1 VERSION
6              
7             This document refers to version PDL::Matrix 0.5 of PDL::Matrix
8              
9             =head1 SYNOPSIS
10              
11             use PDL::Matrix;
12              
13             $m = mpdl [[1,2,3],[4,5,6]];
14             $m = PDL::Matrix->pdl([[1,2,3],[4,5,6]]);
15             $m = msequence(4,3);
16             @dimsa = $x->mdims; # 'dims' is not overloaded
17              
18             $v = vpdl [0,1,2,3]
19             $v = vzeroes(4);
20              
21             =head1 DESCRIPTION
22              
23             =head2 Overview
24              
25             This package tries to help people who want to use PDL for 2D matrix
26             computation with lots of indexing involved. It provides a PDL
27             subclass so one- and two-dimensional piddles that are used as
28             vectors resp and matrices can be typed in using traditional matrix
29             convention.
30              
31             If you want to know more about matrix operation support in PDL, you
32             want to read L or L.
33              
34             The original pdl class refers to the first index as the first row,
35             the second index as the first column of a matrix. Consider
36              
37             print $B = sequence(3,2)
38             [
39             [0 1 2]
40             [3 4 5]
41             ]
42              
43             which gives a 2x3 matrix in terms of the matrix convention, but the
44             constructor used (3,2). This might get more confusing when using
45             slices like sequence(3,2)->slice("1:2,(0)") : with traditional
46             matrix convention one would expect [2 4] instead of [1 2].
47              
48             This subclass PDL::Matrix overloads the constructors and indexing
49             functions of pdls so that they are compatible with the usual matrix
50             convention, where the first dimension refers to the row of a
51             matrix. So now, the above example would be written as
52              
53             print $B = PDL::Matrix->sequence(3,2) # or $B = msequence(3,2)
54             [
55             [0 1]
56             [2 3]
57             [4 5]
58             ]
59              
60             Routines like L or
61             L can be used without any changes.
62              
63             Furthermore one can construct and use vectors as n x 1 matrices
64             without mentioning the second index '1'.
65              
66             =head2 Implementation
67              
68             C works by overloading a number of PDL constructors
69             and methods such that first and second args (corresponding to
70             first and second dims of corresponding matrices) are effectively swapped.
71             It is not yet clear if PDL::Matrix achieves a consistent column-major
72             look-and-feel in this way.
73              
74             =head1 NOTES
75              
76             As of version 0.5 (rewrite by CED) the matrices are stored in the usual
77             way, just constructed and stringified differently. That way indexing
78             and everything else works the way you think it should.
79              
80             =head1 FUNCTIONS
81              
82             =cut
83              
84             package PDL::Matrix;
85              
86             @EXPORT_OK = ();
87              
88              
89             #use PDL::Core;
90             #use PDL::Slatec;
91 1     1   72377 use PDL::Exporter;
  1         2  
  1         7  
92 1     1   6 use Carp;
  1         1  
  1         176  
93              
94             @ISA = qw/PDL::Exporter PDL/;
95              
96             our $VERSION = "0.5";
97             $VERSION = eval $VERSION;
98              
99             #######################################################################=
100             #########
101             #
102             # overloads
103              
104             use overload( '""' => \&string,
105 0     0   0 'x' => sub {my $foo = $_[0]->null();
106 0         0 &PDL::Primitive::matmult(@_[1,0],$foo);
107 0         0 $foo;}
108 1     1   8 );
  1         1  
  1         21  
109              
110             sub string {
111 1     1 0 191 my ($me,@a) = shift;
112 1 50       13 return $me->SUPER::string(@a) unless($me->ndims > 0);
113 0 0       0 $me = $me->dummy(1,1) unless($me->ndims > 1);
114 0         0 $me->xchg(0,1)->SUPER::string(@a);
115             }
116              
117              
118             # --------> constructors
119              
120             =head2 mpdl, PDL::Matrix::pdl
121              
122             =for ref
123              
124             constructs an object of class PDL::Matrix which is a piddle child class.
125              
126             =for example
127              
128             $m = mpdl [[1,2,3],[4,5,6]];
129             $m = PDL::Matrix->pdl([[1,2,3],[4,5,6]]);
130              
131             =cut
132              
133             sub pdl {
134 1     1 1 4 my $class = shift;
135 1         10 my $pdl = $class->SUPER::pdl(@_);
136 1 50       13 if($pdl->ndims > 0) {
137 1 50       5 $pdl = $pdl->dummy(1,1) unless $pdl->ndims > 1;
138 1         66 $pdl = $pdl->xchg(0,1);
139             }
140 1   33     15 bless $pdl, ref $class || $class;
141             }
142              
143             =head2 mzeroes, mones, msequence
144              
145             =for ref
146              
147             constructs a PDL::Matrix object similar to the piddle constructors
148             zeroes, ones, sequence.
149              
150             =cut
151              
152             for my $func (qw /pdl zeroes ones sequence dims/) {
153             push @EXPORT_OK, "m$func";
154 0     0 0 0 eval " sub m$func { PDL::Matrix->$func(\@_) }; ";
  0     0 1 0  
  1     1 1 157  
  0     0 1    
  0     0 1    
155             }
156              
157             =head2 vpdl
158              
159             =for ref
160              
161             constructs an object of class PDL::Matrix which is of matrix
162             dimensions (n x 1)
163              
164             =for example
165              
166             print $v = vpdl [0,1];
167             [
168             [0]
169             [1]
170             ]
171              
172             =cut
173              
174             sub vpdl {
175 0     0 1   my $pdl = PDL->pdl(@_);
176 0           bless $pdl, PDL::Matrix;
177             }
178             push @EXPORT_OK, "vpdl";
179              
180             =head2 vzeroes, vones, vsequence
181              
182             =for ref
183              
184             constructs a PDL::Matrix object with matrix dimensions (n x 1),
185             therefore only the first scalar argument is used.
186              
187             =for example
188              
189             print $v = vsequence(2);
190             [
191             [0]
192             [1]
193             ]
194              
195             =cut
196              
197             for my $func (qw /zeroes ones sequence/) {
198             push @EXPORT_OK, "v$func";
199             my $code = << "EOE";
200              
201             sub v$func {
202             my \@arg = \@_;
203             ref(\$arg[0]) ne 'PDL::Type' ? (\@arg = (\$arg[0],1)) :
204             (\@arg = (\$arg[0],\$arg[1],1));
205             PDL::Matrix->$func(\@arg);
206             }
207              
208             EOE
209             # print "evaluating $code\n";
210 0 0   0 1   eval $code;
  0 0   0 1    
  0 0   0 1    
  0            
  0            
  0            
  0            
  0            
  0            
211             }
212              
213              
214              
215 1     1   205 eval "use PDL::Slatec";
  0         0  
  0         0  
216              
217             my $has_slatec = ($@ ? 0 : 1);
218             sub inv {
219 0     0 0   my $self = shift;
220 0 0         croak "inv: PDL::Slatec not available" unless $has_slatec;
221 0           return $self->matinv;
222             }
223              
224             =head2 kroneckerproduct
225              
226             =for ref
227              
228             returns kroneckerproduct of two matrices. This is not efficiently
229             implemented.
230              
231             =for example
232             print kroneckerproduct(msequence(2,2),mones(2,2))
233             [
234             [0 0 1 1]
235             [0 0 1 1]
236             [2 2 3 3]
237             [2 2 3 3]
238             ]
239              
240             =cut
241              
242             # returns kroneckerproduct of two matrices
243             sub kroneckerproduct {
244 0     0 1   my @arg = @_;
245            
246 0           my ($r0,$c0) = $arg[0]->mdims;
247 0           my ($r1,$c1) = $arg[1]->mdims;
248            
249 0           my $out = mzeroes($r0*$r1,$c0*$c1);
250            
251 0           for (my $i=0;$i<$r0;$i++) {
252 0           for (my $j=0;$j<$c0;$j++) {
253 0           ($_ = $out->slice(($i*$r1).":".(($i+1)*$r1-1).",".
254             ($j*$c1).":".(($j+1)*$c1-1)) ) .= $arg[0]->at($i,$j) * $arg[1];
255             }
256             }
257            
258 0           return $out;
259             }
260             push @EXPORT_OK, "kroneckerproduct";
261              
262             sub rotate {
263 0     0 0   my ($self,@args) = @_;
264 0           return $self->transpose->SUPER::rotate(@args)->transpose;
265             }
266              
267              
268             sub msumover {
269 0     0 0   my ($mpdl) = @_;
270 0           return PDL::sumover(transpose($mpdl)->xchg(0,2));
271             }
272             push @EXPORT_OK, "msumover";
273              
274              
275             =head2 det_general
276              
277             =for ref
278              
279             returns a generalized determinant of a matrix. If the matrix is not
280             regular, one can specify the rank of the matrix and the corresponding
281             subdeterminant is returned. This is implemented using the C
282             function.
283              
284             =for example
285             print msequence(3,3)->determinant(2) # determinant of
286             # regular 2x2 submatrix
287             -24
288              
289             =cut
290              
291             #
292             sub det_general {
293 0     0 1   my ($mpdl,$rank) = @_;
294 0           my $eigenvalues = (PDL::Math::eigens($mpdl))[1];
295 0           my @sort = list(PDL::Ufunc::qsorti(abs($eigenvalues)));
296 0           $eigenvalues = $eigenvalues->dice([@sort[-$rank..-1]]);
297 0           PDL::Ufunc::dprod($eigenvalues);
298             }
299              
300             =head2 trace
301              
302             =for ref
303              
304             returns the trace of a matrix (sum of diagonals)
305              
306             =cut
307              
308             sub trace {
309 0     0 1   my ($mpdl) = @_;
310 0           $mpdl->diagonal(0,1)->sum;
311             }
312              
313             # this has to be overloaded so that the PDL::slice
314             # is called and not PDL::Matrix::slice :-(
315             sub dummy($$;$) {
316 0     0 0   my ($pdl,$dim) = @_;
317 0 0         $dim = $pdl->getndims+1+$dim if $dim < 0;
318 0 0 0       barf ("too high/low dimension in call to dummy, allowed min/max=0/"
319             . $_[0]->getndims)
320             if $dim>$pdl->getndims || $dim < 0;
321 0 0         $_[2] = 1 if ($#_ < 2);
322 0           $pdl->PDL::slice((','x$dim)."*$_[2]");
323             }
324              
325              
326             # now some of my very own helper functions...
327             # stupid function to print a PDL::Matrix object in Maple code
328             sub stringifymaple {
329 0     0 0   my ($self,@args) = @_;
330              
331 0           my ($dimR,$dimC) = mdims($self);
332 0           my $s;
333              
334 0 0         $s .= $args[0].":=" unless $args[0] eq "";
335 0 0         if (defined($dimR)) {
336 0           $s .= "matrix($dimR,$dimC,[";
337 0           for(my $i=0;$i<$dimR;++$i) {
338 0           $s .= "[";
339 0           for(my $j=0;$j<$dimC;++$j) {
340 0           $s .= $self->at($i,$j);
341 0 0         $s .= "," if $j+1<$dimC;
342             }
343 0           $s .= "]";
344 0 0         $s .= "," if $i+1<$dimR;
345             }
346 0           $s .= "])";
347             }
348             else {
349 0           $s = "vector($dimC,[";
350 0           for(my $i=0;$i<$dimC;++$i) {
351 0           $s .= $self->at($i);
352 0 0         $s .= "," if $i+1<$dimC;
353             }
354 0           $s .= "])";
355             }
356 0           return $s;
357             }
358             sub printmaple {
359 0     0 0   print stringifymaple(@_).";\n";
360             }
361              
362             # stupid function to print a PDL::Matrix object in (La)TeX code
363             sub stringifyTeX {
364 0     0 0   my ($self,@args) = @_;
365              
366 0           my ($dimR,$dimC) = mdims($self);
367 0           my $s;
368              
369 0 0         $s .= $args[0]."=" unless $args[0] eq "";
370 0           $s .= "\\begin{pmatrix}\n";
371 0           for(my $i=0;$i<$dimR;++$i) {
372 0           for(my $j=0;$j<$dimC;++$j) {
373 0           $s .= $self->at($i,$j);
374 0 0         $s .= " & " if $j+1<$dimC;
375             }
376 0 0         $s .= " \\\\ \n" if $i+1<$dimR;
377             }
378 0           $s .= "\n \\end{pmatrix}\n";
379              
380 0           return $s;
381             }
382              
383             sub printTeX {
384 0     0 0   print stringifyTeX(@_)."\n";
385             }
386              
387             =pod
388              
389             =begin comment
390              
391             DAL commented this out 17-June-2008. It didn't work, it used the
392             outmoded (and incorrect) ~-is-transpose convention, and it wasn't
393             necessary since the regular cross product worked fine.
394              
395             =head2 vcrossp, PDL::Matrix::crossp
396              
397             =for ref
398              
399             similar to PDL::crossp, however reflecting PDL::Matrix notations
400              
401             #=cut
402              
403             # crossp for my special vectors
404             sub crossp {
405             my ($pdl1,$pdl2) = @_;
406             return PDL::transpose(PDL::crossp(~$pdl1,~$pdl2));
407             }
408             sub vcrossp { PDL::Matrix->crossp(\@_) }
409             push @EXPORT_OK, "vcrossp";
410              
411             =end comment
412              
413             =cut
414              
415             %EXPORT_TAGS = (Func=>[@EXPORT_OK]);
416              
417             1;
418              
419             =head1 BUGS AND PROBLEMS
420              
421             Because we change the way piddles are constructed, not all pdl
422             operators may be applied to piddle-matrices. The inner product is not
423             redefined. We might have missed some functions/methods. Internal
424             consistency of our approach needs yet to be established.
425              
426             Because PDL::Matrix changes the way slicing behaves, it breaks many
427             operators, notably those in MatrixOps.
428              
429             =head1 TODO
430              
431             check all PDL functions, benchmarks, optimization, lots of other things ...
432              
433             =head1 AUTHOR(S)
434              
435             Stephan Heuel (stephan@heuel.org), Christian Soeller
436             (c.soeller@auckland.ac.nz).
437              
438             =head1 COPYRIGHT
439              
440             All rights reserved. There is no warranty. You are allowed to
441             redistribute this software / documentation under certain
442             conditions. For details, see the file COPYING in the PDL
443             distribution. If this file is separated from the PDL distribution, the
444             copyright notice should be included in the file.
445              
446             =cut