152 lines
3.5 KiB
Perl
152 lines
3.5 KiB
Perl
use feature 'say';
|
|
use ntheory qw/forsetproduct/;
|
|
|
|
# Deconvolution of N dimensional matrices
|
|
sub deconvolve_N {
|
|
our @g; local *g = shift;
|
|
our @f; local *f = shift;
|
|
my @df = shape(@f);
|
|
my @dg = shape(@g);
|
|
my @hsize;
|
|
push @hsize, $dg[$_] - $df[$_] + 1 for 0..$#df;
|
|
my @toSolve = map { [row(\@g, \@f, \@hsize, $_)] } coords(shape(@g));
|
|
rref( \@toSolve );
|
|
|
|
my @h;
|
|
my $n = 0;
|
|
for (coords(@hsize)) {
|
|
my($k,$j,$i) = split ' ', $_;
|
|
$h[$i][$j][$k] = $toSolve[$n++][-1];
|
|
}
|
|
@h;
|
|
}
|
|
|
|
sub row {
|
|
our @g; local *g = shift;
|
|
our @f; local *f = shift;
|
|
our @hsize; local *hsize = shift;
|
|
my @gc = reverse split ' ', shift;
|
|
|
|
my @row;
|
|
my @fdim = shape(@f);
|
|
for (coords(@hsize)) {
|
|
my @hc = reverse split ' ', $_;
|
|
my @fc;
|
|
for my $i (0..$#hc) {
|
|
my $window = $gc[$i] - $hc[$i];
|
|
push(@fc, $window), next if 0 <= $window && $window < $fdim[$i];
|
|
}
|
|
push @row, $#fc == $#hc ? $f [$fc[0]] [$fc[1]] [$fc[2]] : 0;
|
|
}
|
|
push @row, $g [$gc[0]] [$gc[1]] [$gc[2]];
|
|
return @row;
|
|
}
|
|
|
|
sub rref {
|
|
our @m; local *m = shift;
|
|
@m or return;
|
|
my ($lead, $rows, $cols) = (0, scalar(@m), scalar(@{$m[0]}));
|
|
|
|
foreach my $r (0 .. $rows - 1) {
|
|
$lead < $cols or return;
|
|
my $i = $r;
|
|
|
|
until ($m[$i][$lead])
|
|
{++$i == $rows or next;
|
|
$i = $r;
|
|
++$lead == $cols and return;}
|
|
|
|
@m[$i, $r] = @m[$r, $i];
|
|
my $lv = $m[$r][$lead];
|
|
$_ /= $lv foreach @{ $m[$r] };
|
|
|
|
my @mr = @{ $m[$r] };
|
|
foreach my $i (0 .. $rows - 1)
|
|
{$i == $r and next;
|
|
($lv, my $n) = ($m[$i][$lead], -1);
|
|
$_ -= $lv * $mr[++$n] foreach @{ $m[$i] };}
|
|
|
|
++$lead;}
|
|
}
|
|
|
|
# Constructs an AoA of coordinates to all elements of N dimensional array
|
|
sub coords {
|
|
my(@dimensions) = reverse @_;
|
|
my(@ranges,@coords);
|
|
push @ranges, [0..$_-1] for @dimensions;
|
|
forsetproduct { push @coords, "@_" } @ranges;
|
|
@coords;
|
|
}
|
|
|
|
sub shape {
|
|
my(@dim);
|
|
push @dim, scalar @_;
|
|
push @dim, shape(@{$_[0]}) if 'ARRAY' eq ref $_[0];
|
|
@dim;
|
|
}
|
|
|
|
# Pretty printer for N dimensional arrays
|
|
# Assumes if first element in level is an array, then all are
|
|
sub pretty_print {
|
|
my($i, @a) = @_;
|
|
if ('ARRAY' eq ref $a[0]) {
|
|
say ' 'x$i, '[';
|
|
pretty_print($i+2, @$_) for @a;
|
|
say ' 'x$i, ']', $i ? ',' : '';
|
|
} else {
|
|
say ' 'x$i, '[', sprintf("@{['%5s'x@a]}",@a), ']', $i ? ',' : '';
|
|
}
|
|
}
|
|
|
|
my @f = (
|
|
[
|
|
[ -9, 5, -8 ],
|
|
[ 3, 5, 1 ],
|
|
],
|
|
[
|
|
[ -1, -7, 2 ],
|
|
[ -5, -6, 6 ],
|
|
],
|
|
[
|
|
[ 8, 5, 8 ],
|
|
[ -2, -6, -4 ],
|
|
]
|
|
);
|
|
|
|
my @g = (
|
|
[
|
|
[ 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 ],
|
|
]
|
|
);
|
|
|
|
my @h = deconvolve_N( \@g, \@f );
|
|
my @ff = deconvolve_N( \@g, \@h );
|
|
|
|
my $d = scalar shape(@g);
|
|
print "${d}D arrays:\n";
|
|
print "h =\n";
|
|
pretty_print(0,@h);
|
|
print "\nff =\n";
|
|
pretty_print(0,@ff);
|