140 lines
3.5 KiB
Raku
140 lines
3.5 KiB
Raku
# 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 );
|