2D
Translation of Tcl.
# Deconvolution of N dimensional matrices.
sub deconvolve-N ( @g, @f ) {
my @hsize = @g.shape »-« @f.shape »+» 1;
my @toSolve = coords(@g.shape).map:
{ [row(@g, @f, $^coords, @hsize)] };
my @solved = rref( @toSolve );
my @h;
for flat coords(@hsize) Z @solved[*;*-1] -> $_, $v {
@h.AT-POS(|$_) = $v;
}
@h
}
# Construct a row for each value in @g to be sent to the simultaneous equation solver
sub row ( @g, @f, @gcoord, $hsize ) {
my @row;
@gcoord = @gcoord[(^@f.shape)]; # clip extraneous values
for coords( $hsize ) -> @hc {
my @fcoord;
for ^@hc -> $i {
my $window = @gcoord[$i] - @hc[$i];
@fcoord.push($window) and next if 0 ≤ $window < @f.shape[$i];
last;
}
@row.push: @fcoord == @hc ?? @f.AT-POS(|@fcoord) !! 0;
}
@row.push: @g.AT-POS(|@gcoord);
@row
}
# Constructs an AoA of coordinates to all elements of N dimensional array
sub coords ( @dim ) {
@[reverse $_ for [X] ([^$_] for reverse @dim)]
}
# Reduced Row Echelon Form simultaneous equation solver
# Can handle over-specified systems (N unknowns in N + M equations)
sub rref (@m) {
@m = trim-system @m;
my ($lead, $rows, $cols) = 0, @m, @m[0];
for ^$rows -> $r {
return @m unless $lead < $cols;
my $i = $r;
until @m[$i;$lead] {
next unless ++$i == $rows;
$i = $r;
return @m if ++$lead == $cols;
}
@m[$i, $r] = @m[$r, $i] if $r != $i;
@m[$r] »/=» $ = @m[$r;$lead];
for ^$rows -> $n {
next if $n == $r;
@m[$n] »-=» @m[$r] »×» (@m[$n;$lead] // 0);
}
++$lead;
}
@m
}
# Reduce to N equations in N unknowns; a no-op unless rows > cols
sub trim-system (@m) {
return @m unless @m ≥ @m[0];
my (\vars, @t) = @m[0] - 1;
for ^vars -> \lead {
for ^@m -> \row {
@t.append: @m.splice(row, 1) and last if @m[row;lead];
}
}
while @t < vars and @m { @t.push: shift @m }
@t
}
# Pretty printer for N dimensional arrays
# Assumes if first element in level is an array, then all are
sub pretty-print ( @array, $indent = 0 ) {
if @array[0] ~~ Array {
say ' ' x $indent,"[";
pretty-print( $_, $indent + 2 ) for @array;
say ' ' x $indent, "]{$indent??','!!''}";
} else {
say ' ' x $indent, "[{say_it(@array)} ]{$indent??','!!''}";
}
sub say_it ( @array ) { return join ",", @array».fmt("%4s"); }
}
my @f[3;2;3] = (
[
[ -9, 5, -8 ],
[ 3, 5, 1 ],
],
[
[ -1, -7, 2 ],
[ -5, -6, 6 ],
],
[
[ 8, 5, 8 ],
[ -2, -6, -4 ],
]
);
my @g[4;4;6] = (
[
[ 54, 42, 53, -42, 85, -72 ],
[ 45,-170, 94, -36, 48, 73 ],
[ -39, 65,-112, -16, -78, -72 ],
[ 6, -11, -6, 62, 49, 8 ],
],
[
[ -57, 49, -23, 52,-135, 66 ],
[ -23, 127, -58, -5,-118, 64 ],
[ 87, -16, 121, 23, -41, -12 ],
[ -19, 29, 35,-148, -11, 45 ],
],
[
[ -55,-147,-146, -31, 55, 60 ],
[ -88, -45, -28, 46, -26,-144 ],
[ -12,-107, -34, 150, 249, 66 ],
[ 11, -15, -34, 27, -78, -50 ],
],
[
[ 56, 67, 108, 4, 2, -48 ],
[ 58, 67, 89, 32, 32, -8 ],
[ -42, -31,-103, -30, -23, -8 ],
[ 6, 4, -26, -10, 26, 12 ],
]
);
say "# {+@f.shape}D array:";
my @h = deconvolve-N( @g, @f );
say "h =";
pretty-print( @h );
my @h-shaped[2;3;4] = @(deconvolve-N( @g, @f ));
my @ff = deconvolve-N( @g, @h-shaped );
say "\nff =";
pretty-print( @ff );
Output:
# 3D array:
h =
[
[
[ -6, -8, -5, 9 ],
[ -7, 9, -6, -8 ],
[ 2, -7, 9, 8 ],
],
[
[ 7, 4, 4, -6 ],
[ 9, 9, 4, -4 ],
[ -3, 7, -2, -3 ],
],
]
ff =
[
[
[ -9, 5, -8 ],
[ 3, 5, 1 ],
],
[
[ -1, -7, 2 ],
[ -5, -6, 6 ],
],
[
[ 8, 5, 8 ],
[ -2, -6, -4 ],
],
]
Last updated