transpose degree n matrices

Author: Alex-Jordan

lib/Value/Matrix.pm

@@ -848,29 +848,119 @@ sub twiddle {
 	return $self->make(@coords);
 }
 
+=head3 C<slice>
+
+Produce the degree (n-1) Matrix defined by a given index and value for that index. If n is 1,
+this produces a Real/Complex/Fraction.
+
+Usage:
+
+    $A = Matrix([ [ 1, 2, 3, 4 ], [ 5, 6, 7, 8 ], [ 9, 10, 11, 12 ] ]);
+    $A->slice(1, 2)            # will be same as Matrix([5, 6, 7, 8])
+    $A->slice(2, 3)            # will be same as Matrix([3, 7, 11])
+
+    $B = Matrix([ [ [ 1, 2 ], [ 3, 4 ] ], [ [ 5, 6 ], [ 7, 8 ] ] ]);
+    $B->slice(1, 2)            # will be same as Matrix([ [ 5, 6 ], [ 7, 8 ] ])
+    $B->slice(2, 1)            # will be same as Matrix([ [ 1, 2 ], [ 5, 6 ] ])
+    $B->slice(3, 1)            # will be same as Matrix([ [ 1, 3 ], [ 5, 7 ] ])
+
+=cut
+
+sub slice {
+	my $self = shift;
+	my ($index, $value) = @_;
+	my @d = $self->dimensions;
+	my $d = scalar(@d);
+	my $w = $d[0];
+	Value::Error("index must be an integer from 1 to $d") unless ($index == int($index) && $index >= 1 && $index <= $d);
+	my $M = $self->data;
+	if ($index == 1) {
+		Value::Error("value must be an integer from 1 to $w")
+			unless ($value == int($value) && $value >= 1 && $value <= $w);
+		return $M->[ $value - 1 ];
+		return $self->make($M->[ $value - 1 ]);
+	} else {
+		my @rows;
+		for (1 .. $w) {
+			push @rows, $M->[ $_ - 1 ]->slice($index - 1, $value);
+		}
+		return $self->make(@rows);
+	}
+}
+
 =head3 C<transpose>
 
-Take the transpose of a matrix.
+Take the transpose of a matrix. For a degree 1 Matrix, first promote to a degree 2 Matrix.
+For a degree n Matrix, apply a permutation of the indices. The default permutation transposes the
+last two indices. To specify a permutation, provide an array reference representing a cycle
+or an array of array references that represents a product of cycles. If a permutation is not
+specified, the default is the usual transposition of the last two indices.
 
 Usage:
 
     $A = Matrix([ [ 1, 2, 3, 4 ], [ 5, 6, 7, 8 ], [ 9, 10, 11, 12 ] ]);
-    $A->transpose;
+    $A->transpose              # will be [ [ 1, 5, 9 ], [ 2, 6, 10 ], [ 3, 7, 11 ], [ 4, 8, 12 ] ]
 
+    $B = Matrix([ [ [ 1, 2 ], [ 3, 4 ] ], [ [ 5, 6 ], [ 7, 8 ] ] ]);
+    $B->transpose([1, 2, 3])   # will be [ [ [ 1, 3 ], [ 5, 7 ] ], [ [2 , 4 ], [ 6, 8 ] ] ]
+
+    $C = Matrix([ [ [ [ 1, 2 ], [ 3, 4 ] ], [ [ 5, 6 ], [ 7, 8 ] ] ], [ [ [ 9, A ], [ B, C ] ], [ [ D, E ], [ F, 0 ] ] ] ]);
+    $C->transpose([ [ 1, 2], [3, 4] ])   # will be [ [ [ [ 1, 3 ], [ 2, 4 ] ], [ [ 9, B ], [ A, C ] ] ], [ [ [ 5, 7 ], [ 6, 8 ] ], [ [ D, F ], [ E, 0 ] ] ]
 =cut
 
 sub transpose {
-	my $self = promote(@_);
+	my $self = shift;
+	my $p    = shift;
 	my @d    = $self->dimensions;
-	if (scalar(@d) == 1) { @d = (1, @d); $self = $self->make($self) }
-	Value::Error("Can't transpose %d-dimensional matrices", scalar(@d)) unless scalar(@d) == 2;
+	my $N    = scalar(@d);
+
+	# elevate a degree 1 Matrix to degree 2
+	if ($N == 1) { @d = (1, @d); $N = 2; $self = $self->make($self) }
+
+	# default to transpose last two indices
+	$p = [ [ $N - 1, $N ] ] unless $p;
+
+	# build the permutation hash from cycles
+	my %p;
+	if (ref $p eq 'HASH') {
+		%p = %{$p};
+	} else {
+		$p = [$p] unless ref($p->[0]);
+		my @p = (1 .. $N);
+		for my $cycle (@{$p}) {
+			next unless defined $cycle->[0];
+			my $tmp = $p[ $cycle->[0] - 1 ];
+			for my $i (0 .. $#{$cycle} - 1) {
+				$p[ $cycle->[$i] - 1 ] = $p[ $cycle->[ $i + 1 ] - 1 ];
+			}
+			$p[ $cycle->[ $#{$cycle} ] - 1 ] = $tmp;
+		}
+		%p = map { $_ => $p[ $_ - 1 ] } (1 .. $N);
+	}
+	%p = reverse %p;
 
 	my @M = ();
-	my $M = $self->data;
-	for my $j (0 .. $d[1] - 1) {
-		my @row = ();
-		for my $i (0 .. $d[0] - 1) { push(@row, $M->[$i]->data->[$j]) }
-		push(@M, $self->make(@row));
+	if ($N == 2) {
+		return $self if ($p{1} == 1);
+		my $M = $self->data;
+		for my $j (0 .. $d[1] - 1) {
+			my @row = ();
+			for my $i (0 .. $d[0] - 1) { push(@row, $M->[$i]->data->[$j]) }
+			push(@M, $self->make(@row));
+		}
+	} else {
+		# reduce the permutation hash
+		my @q  = map { $p{$_} } (1 .. $N);
+		my $p1 = shift @q;
+		for (@q) {
+			$_-- if ($_ >= $p1);
+		}
+		my %q = map { $_ => $q[ $_ - 1 ] } (1 .. $N - 1);
+
+		for my $j (1 .. $d[ $p1 - 1 ]) {
+			my $slice = $self->slice($p1, $j);
+			push(@M, $slice->class eq 'Matrix' ? $slice->transpose(\%q) : $slice);
+		}
 	}
 	return $self->make(@M);
 }