 Matrix.pm.orig Wed Oct 24 09:20:10 2001
+++ Matrix.pm Sun Oct 14 14:48:40 2001
@@ 32,7 +32,7 @@
The original pdl class refers to the first index as the first row,
the second index as the first column of a matrix. Consider
 print $b = sequence(3,2)
+ print $B = sequence(3,2)
[
[0 1 2]
[3 4 5]
@@ 48,7 +48,7 @@
convention, where the first dimension refers to the row of a
matrix. So now, the above example would be written as
 print $b = PDL::Matrix>sequence(3,2)
+ print $B = PDL::Matrix>sequence(3,2) # or $B = msequence(3,2)
[
[0 1]
[2 3]
@@ 71,133 +71,20 @@
=head1 FUNCTIONS
=head2 mpdl, PDL::Matrix::pdl

=for ref

constructs an object of class PDL::Matrix which is a piddle child
class, where the first index refers to the first column of the
twodimensional piddle.

=for example

 $m = mpdl [[1,2,3],[4,5,6]];
 $m = PDL::Matrix>pdl([[1,2,3],[4,5,6]]);

=head2 mzeroes, mones, msequence

=for ref

constructs a PDL::Matrix object similar to the piddle constructors
zeroes, ones, sequence

=head2 vpdl

=for ref

constructs an object of class PDL::Matrix which is of matrix
dimensions (n x 1)

=for example

 print $v = vpdl [0,1];
 [
 [0]
 [1]
 ]

=head2 vzeroes, vones, vsequence

=for ref

constructs a PDL::Matrix object with matrix dimensions (n x 1),
therefore only the first scalar argument is used.


=head2 PDL::Matrix::slice, PDL::Matrix::dice

=for ref

same as slice, dice for normal piddles, but reflecting the matrix
convention by swapping the first two arguments.

=for example

 print sequence(3,2)>slice("1:2,(0)") # piddle
 [1 2]
 print msequence(3,2)>slice("1:2,(0)") # PDL::Matrix
 [2 4]

=head2 PDL::Matrix::at

=for ref

same as at for piddles, but reflecting the matrix convention by
swapping the first two arguments

If only one scalar argument is used, we assume the object to be a
vector and look only at the first column.

=head2 mdims

=for ref

returns the dimensions of the PDL::Matrix object in matrix
convention

C is NOT overloaded by PDL::Matrix to make sure that
methods like PDL::transpose still work. So use C to get
the dims in the PDL::Matrix notation.

=for example

 print msequence(3,2)>mdims
 3 2

=head2 vcrossp, PDL::Matrix::crossp

=for ref

similar to PDL::crossp, however reflecting PDL::Matrix notations

=head1 BUGS AND PROBLEMS

Because we change the way piddles are constructed, not all pdl
operators may be applied to piddlematrices. The inner product is not
redefined. We might have missed some functions/methods. Internal
consistency of our approach needs yet to be established.

=head1 TODO

benchmarks, optimization ...

=head1 AUTHOR(S)

Stephan Heuel (stephan@heuel.org) with lots of help from Christian
Soeller (c.soeller@auckland.ac.nz).

=head1 COPYRIGHT

All rights reserved. There is no warranty. You are allowed to
redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL
distribution. If this file is separated from the PDL distribution, the
copyright notice should be included in the file.

=cut

package PDL::Matrix;
@EXPORT_OK = ();
use PDL::Core;
+
+#use PDL::Core;
use PDL::Slatec;
use PDL::Exporter;
@ISA = qw/PDL::Exporter PDL/;
$VERSION = "0.01";
+$VERSION = "0.02";
#######################################################################=
#########
@@ 206,13 +93,36 @@
# > constructors
# pdl is basically the same
+=head2 mpdl, PDL::Matrix::pdl
+
+=for ref
+
+constructs an object of class PDL::Matrix which is a piddle child
+class, where the first index refers to the first column of the
+twodimensional piddle.
+
+=for example
+
+ $m = mpdl [[1,2,3],[4,5,6]];
+ $m = PDL::Matrix>pdl([[1,2,3],[4,5,6]]);
+
+=cut
+
sub pdl {
my $class = shift;
my $pdl = $class>SUPER::pdl(@_);
bless $pdl, ref $class  $class;
}
+=head2 mzeroes, mones, msequence
+
+=for ref
+
+constructs a PDL::Matrix object similar to the piddle constructors
+zeroes, ones, sequence
+
+=cut
+
# for constructors with specified dimensions I only have to swap these
# dimensions
for my $func (qw /zeroes ones sequence/) {
@@ 248,8 +158,80 @@
eval $code;
}
# > methods
# the slice arguments have to be changed to reflect the swapping
+=head2 vpdl
+
+=for ref
+
+constructs an object of class PDL::Matrix which is of matrix
+dimensions (n x 1)
+
+=for example
+
+ print $v = vpdl [0,1];
+ [
+ [0]
+ [1]
+ ]
+
+=cut
+
+sub vpdl {
+ my $pdl = transpose(PDL>pdl(@_));
+ bless $pdl, PDL::Matrix;
+}
+push @EXPORT_OK, "vpdl";
+
+=head2 vzeroes, vones, vsequence
+
+=for ref
+
+constructs a PDL::Matrix object with matrix dimensions (n x 1),
+therefore only the first scalar argument is used.
+
+=for example
+
+ print $v = vsequence(2);
+ [
+ [0]
+ [1]
+ ]
+
+=cut
+
+for my $func (qw /zeroes ones sequence/) {
+ push @EXPORT_OK, "v$func";
+ my $code = << "EOE";
+
+sub v$func {
+ my \@arg = \@_;
+ ref(\$arg[0]) ne 'PDL::Type' ? (\@arg = (\$arg[0],1)) :
+ (\@arg = (\$arg[0],\$arg[1],1));
+ PDL::Matrix>$func(\@arg);
+}
+
+EOE
+# print "evaluating $code\n";
+ eval $code;
+}
+
+
+
+=head2 PDL::Matrix::slice, PDL::Matrix::dice
+
+=for ref
+
+same as slice, dice for normal piddles, but reflecting the matrix
+convention by swapping the first two arguments.
+
+=for example
+
+ print sequence(3,2)>slice("1:2,(0)") # piddle
+ [1 2]
+ print msequence(3,2)>slice("1:2,(0)") # PDL::Matrix
+ [2 4]
+
+=cut
+
sub slice {
my $self = shift;
my $ind = shift;
@@ 272,6 +254,17 @@
}
+=head2 PDL::Matrix::at
+
+=for ref
+
+same as at for piddles, but reflecting the matrix convention by
+swapping the first two arguments
+
+If only one scalar argument is used, we assume the object to be a
+vector and look only at the first column.
+
+=cut
# swap arguments if number of arguments is greater than 1
# if its one, look at the first column, assuming it is a vector
@@ 283,6 +276,61 @@
$self>SUPER::at(@arg);
}
+=head2 PDL::Matrix::set
+
+=for ref
+
+set a particular value in a PDL::Matrix object. Note that this has to
+be called as an object method rather than a function
+
+=for example
+
+print msequence(3,3)>set(2,0,1) # ok with PDL::Matrix convention
+[
+ [ 0 1 2]
+ [ 3 4 5]
+ [1 7 8]
+]
+
+print set msequence(3,3), 2,0,1 # does not conform with PDL::Matrix convention
+[
+ [ 0 1 1]
+ [ 3 4 5]
+ [ 6 7 8]
+]
+
+
+=cut
+
+# this does not work as the usual set, but can be used via
+# $M>set(0,2,0.4) which sets the element in the first row and third
+# column to 0.4
+sub set{ # Sets a particular single value
+ barf("Usage: set($pdl, $x, $y,.., $value)") if $#_<2;
+ my $self = shift; my $value = pop @_;
+ my @arg = @_;
+ @arg >=2 ? @arg[1,0] = @arg[0,1] : ( (@arg == 1 && $self>dims == 2) ?
+ @arg = (0,$arg[0]) : 1 );
+ PDL::Core::set_c ($self, [@arg], $value);
+ return $self;
+}
+
+=head2 PDL::Matrix::reshape
+
+=for ref
+
+same as reshape for piddles, but reflecting the matrix convention by
+swapping the first two arguments
+
+=cut
+
+sub reshape {
+ my $self = shift;
+ my @arg = @_;
+ @arg >=2 ? @arg[1,0] = @arg[0,1] : (@arg = (1,$arg[0]));
+ $self>SUPER::reshape(@arg);
+}
+
# this is needed because only doing
# > use overload '~' => \&PDL::transpose;
# would not give an object of type PDL::Matrix back
@@ 295,6 +343,24 @@
use overload '~' => \&PDL::Matrix::transpose;
+=head2 mdims
+
+=for ref
+
+returns the dimensions of the PDL::Matrix object in matrix
+convention
+
+C is NOT overloaded by PDL::Matrix to make sure that
+methods like PDL::transpose still work. So use C to get
+the dims in the PDL::Matrix notation.
+
+=for example
+
+ print msequence(3,2)>mdims
+ 3 2
+
+=cut
+
# I cannot overload dims, because it is apparently used many times
# in methods like transpose! :(
# this is not nice but I don't know how to avoid this
@@ 310,6 +376,96 @@
return matinv($self);
}
+
+=head2 kroneckerproduct
+
+=for ref
+
+returns kroneckerproduct of two matrices. This is not efficiently
+implemented.
+
+=for example
+print kroneckerproduct(msequence(2,2),mones(2,2))
+[
+ [0 0 1 1]
+ [0 0 1 1]
+ [2 2 3 3]
+ [2 2 3 3]
+]
+
+=cut
+
+# returns kroneckerproduct of two matrices
+sub kroneckerproduct {
+ my @arg = @_;
+
+ my ($r0,$c0) = $arg[0]>mdims;
+ my ($r1,$c1) = $arg[1]>mdims;
+
+ my $out = mzeroes($r0*$r1,$c0*$c1);
+
+ for (my $i=0;$i<$r0;$i++) {
+ for (my $j=0;$j<$c0;$j++) {
+ ($_ = $out>slice(($i*$r1).":".(($i+1)*$r11).",".
+ ($j*$c1).":".(($j+1)*$c11)) ) .= $arg[0]>at($i,$j) * $arg[1];
+ }
+ }
+
+ return $out;
+}
+push @EXPORT_OK, "kroneckerproduct";
+
+sub rotate {
+ my ($self,@args) = @_;
+ return $self>transpose>SUPER::rotate(@args)>transpose;
+}
+
+
+sub msumover {
+ my ($mpdl) = @_;
+ return PDL::sumover(transpose($mpdl)>xchg(0,2));
+}
+push @EXPORT_OK, "msumover";
+
+
+=head2 det_general
+
+=for ref
+
+returns a generalized determinant of a matrix. If the matrix is not
+regular, one can specify the rank of the matrix and the corresponding
+subdeterminant is returned. This is implemented using the C
+function.
+
+=for example
+print msequence(3,3)>determinant(2) # determinant of
+ # regular 2x2 submatrix
+24
+
+=cut
+
+#
+sub det_general {
+ my ($mpdl,$rank) = @_;
+ my $eigenvalues = (PDL::Math::eigens($mpdl))[1];
+ my @sort = list(PDL::Ufunc::qsorti(abs($eigenvalues)));
+ $eigenvalues = $eigenvalues>dice([@sort[$rank..1]]);
+ PDL::Ufunc::dprod($eigenvalues);
+}
+
+=head2 trace
+
+=for ref
+
+returns the trace of a matrix (sum of diagonals)
+
+=cut
+
+sub trace {
+ my ($mpdl) = @_;
+ $mpdl>diagonal(0,1)>sum;
+}
+
# this has to be overloaded so that the PDL::slice
# is called and not PDL::Matrix::slice :(
sub dummy($$;$) {
@@ 323,6 +479,7 @@
}
+# now some of my very own helper functions...
# stupid function to print a PDL::Matrix object in Maple code
sub stringifymaple {
my ($self,@args) = @_;
@@ 358,7 +515,7 @@
print stringifymaple(@_).";\n";
}
# stupid function to print a PDL::Matrix object in Maple code
+# stupid function to print a PDL::Matrix object in (La)TeX code
sub stringifyTeX {
my ($self,@args) = @_;
@@ 384,31 +541,13 @@
}
+=head2 vcrossp, PDL::Matrix::crossp
# functions that construct a vector pdl and that can be
# exported,
for my $func (qw /zeroes ones sequence/) {
 push @EXPORT_OK, "v$func";
 my $code = << "EOE";

sub v$func {
 my \@arg = \@_;
 ref(\$arg[0]) ne 'PDL::Type' ? (\@arg = (\$arg[0],1)) :
 (\@arg = (\$arg[0],\$arg[1],1));
 PDL::Matrix>$func(\@arg);
}
+=for ref
EOE
# print "evaluating $code\n";
 eval $code;
}
+similar to PDL::crossp, however reflecting PDL::Matrix notations
# vpdl
sub vpdl {
 my $pdl = transpose(PDL>pdl(@_));
 bless $pdl, PDL::Matrix;
}
push @EXPORT_OK, "vpdl";
+=cut
# crossp for my special vectors
sub crossp {
@@ 421,3 +560,29 @@
%EXPORT_TAGS = (Func=>[@EXPORT_OK]);
1;
+
+=head1 BUGS AND PROBLEMS
+
+Because we change the way piddles are constructed, not all pdl
+operators may be applied to piddlematrices. The inner product is not
+redefined. We might have missed some functions/methods. Internal
+consistency of our approach needs yet to be established.
+
+=head1 TODO
+
+check all PDL functions, benchmarks, optimization, lots of other things ...
+
+=head1 AUTHOR(S)
+
+Stephan Heuel (stephan@heuel.org), Christian Soeller
+(c.soeller@auckland.ac.nz).
+
+=head1 COPYRIGHT
+
+All rights reserved. There is no warranty. You are allowed to
+redistribute this software / documentation under certain
+conditions. For details, see the file COPYING in the PDL
+distribution. If this file is separated from the PDL distribution, the
+copyright notice should be included in the file.
+
+=cut