Deconvolution/2D+: Difference between revisions

→‎{{header|Perl 6}}: modernize with shaped arrays, more idiomatic code
m (→‎{{header|Perl 6}}: Combine files to make runnable code)
(→‎{{header|Perl 6}}: modernize with shaped arrays, more idiomatic code)
Line 631:
 
Translation of Tcl.
<lang perl6># Deconvolution of N dimensional matriciesmatrices.
sub deconv_NDdeconvolve-N ( @g, @f ) {
my @gsizehsize = size_of @g.shape »-« @f.shape »+» 1;
 
my @fsize = size_of @f;
my @hsizetoSolve = coords(@gsize >>-<< @fsize >>+>> 1;g.shape).map:
{ [row(@g, @f, $^coords, @hsize)] };
 
my @toSolve = loopcoords(@gsize).map:
{ [row(@g, @f, @gsize, $^coords, @fsize, @hsize)] };
my @solved = rref( @toSolve );
 
# Uncomment if you want to see the rref system of equations.
# pretty_print( @solved );
my @h;
for flat coords(@hsize) Z @solved[*;*-1] -> $_, $v {
my $index = 0;
@h.AT-POS(|$_) = $v;
insert(@h, $_, @solved[$index++][*-1]) for loopcoords(@hsize);
return @h;
# Inserts a value in the correct spot in an N dimensional array.
sub insert ( $array is rw, @coords is copy, $value ) {
my $level = @coords.shift;
if +@coords {
insert( $array[$level], @coords, $value );
} else {
$array[$level] = $value;
}
}
return @h;
}
 
# Construct a row for each value in @g to be sent to the simultaneous equation solver
# Returns a list containing the number of elements in
sub row ( @g, @f, @gcoord, $hsize ) {
# each level of an N dimensional array.
sub size_of ( $m is copy ) {
my @size;
while $m ~~ Array { @size.push(+$m); $m = $m[0]; }
return @size;
}
 
# Construct a row (equation) for each value in @g to be sent
# to the simultaneous equation solver.
# @Xsize = Dimensions of @X, # of elems per level.
# @Xcoords = Path to each element of @X given as a series of indicies.
sub row ( @g, @f, @gsize, @gcoords, @fsize, $hsize ) {
my @row;
@gcoord = @gcoord[(^@f.shape)]; # clip extraneous values
for loopcoords( $hsize ) -> @hcoords {
for coords( $hsize ) my-> @fcoords;hc {
formy ^@hcoords -> $index {fcoord;
for my $window = ^@gcoords[$index]hc -> @hcoords[$index];i {
@fcoords.push(my $window) and next if 0 <= @gcoord[$windowi] <- @fsizehc[$indexi];
@fcoord.push($window) and next if 0 <= $window < @f.shape[$i];
last;
last;
}
@row.push(: +@fcoordsfcoord == +@hcoordshc ?? fetch( @f, .AT-POS(|@fcoords fcoord) !! 0 );
}
@row.push( fetch(: @g, .AT-POS(|@gcoords ) gcoord);
return @row;
 
# Returns the value found in @array with the
# coordinates given in the list of @indicies.
sub fetch (@array, *@indicies) {
my $index = @indicies.shift;
return @array[*-1] ~~ Array
?? fetch( @array[$index], @indicies )
!! @array[$index];
}
}
 
# Constructs an arrayAoA of arrayscoordinates to all elements of coordinatesN todimensional eacharray
sub coords ( @dim ) {
# element in an N dimensional array.
@[reverse $_ for [X] ([^$_] for reverse @dim)];
sub loopcoords ( @hsize ) {
my @hcoords;
for ^([*] @hsize) -> $index {
my @coords;
my $j = $index;
for @hsize -> $dim {
@coords.push( $j % $dim );
$j div= $dim;
}
@hcoords.push( [@coords] );
}
return @hcoords;
}
 
# Reduced Row Echelon Form simultaneous equation solver.
# Can handle over-specified systems of(N unknowns in N + M equations.)
sub rref (@m) {
# (n unknowns in n + m equations)
return unless @m;
sub rref ($m is rw) {
my ($lead, $rows, $cols) = 0, +@m, +@m[0];
return unless $m;
my ($lead, $rows, $cols) = 0, +$m, +$m[0];
 
# Trim off over specified rows if they exist., for efficiency
# Not strictly necessary, but can save a lot of
# redundant calculations.
if $rows >= $cols {
$@m = trim_system($@m);
$rows = +$@m;
}
 
for ^$rows -> $r {
$lead < $cols or return $@m;
my $i = $r;
until $@m[$i][;$lead] {
++$i == $rows or next;
$i = $r;
++$lead == $cols and return $@m;
}
$@m[$i, $r] = $@m[$r, $i] if $r != $i;
my $lv = $@m[$r][;$lead];
$@m[$r] >>»/=>>» $lv;
for ^$rows -> $n {
next if $n == $r;
$@m[$n] >>»-=>>» $@m[$r] >>»*>>» $(@m[$n][;$lead] // 0);
}
++$lead;
}
return $@m;
 
# Reduce a system of equations to nN equations with nN unknowns.
sub trim_system ($m) {
# Looks for an equation with a true value for each position.
my ($vars, @t) = +$m[0]-1;
# If it can't find one, assumes that it has already taken one
# and pushes in the first equation it sees. This assumtion
# will alway be successful except in some cases where an
# under-specified system has been supplied, in which case,
# it would not have been able to reduce the system anyway.
sub trim_system ($m is rw) {
my ($vars, @t) = +$m[0]-1, ();
for ^$vars -> $lead {
for ^$m -> $row {
@t.push(: | $m.splice( $row, 1 ) ) and last if $m[$row][;$lead];
}
}
while (+@t < $vars) and +$m { @t.push(: $m.splice( 0, 1 ) ) };
return @t;
}
}
 
# Pretty printer for N dimensional arrays
# Use with a pretty printer as follows:
# Assumes if first element in level is an array, then all are
 
# Pretty printer for N dimensional arrays. Assumes that
# if the FIRST element in any particular level is an array,
# then ALL the elements at that level are arrays.
sub pretty_print ( @array, $indent = 0 ) {
my $tab = 2;
if @array[0] ~~ Array {
say ' ' x $indent,"[";
pretty_print( $_, $indent + $tab2 ) for @array;
say ' ' x $indent, "]{$indent??','!!''}";
} else {
Line 781 ⟶ 724:
}
 
sub say_it ( @array ) { return join ",", @array».fmt("%4s"); }
return join ",", @array>>.fmt("%4s");
}
}
 
my @f[3;2;3] = (
[
[ -9, 5, -8 ], [ 3, 5, 1 ],
[ 3, 5, 1 ],
],
[
[ -1, -7, 2 ], [ -5, -6, 6 ],
[ -5, -6, 6 ],
],
[
[ 8, 5, 8 ], [ -2, -6, -4 ],
[ -2, -6, -4 ],
]
);
 
my @g[4;4;6] = (
[
[ 54, 42, 53, -42, 85, -72 ],
Line 825 ⟶ 769:
);
 
say "# {+@f.shape}D array:";
=begin skip_output
 
say "g =";
pretty_print( @g );
 
say '-' x 79;
 
say "f =";
pretty_print( @f );
 
say '-' x 79;
=end skip_output
 
say "# {+size_of(@f)}D array:";
say "h =";
pretty_print( deconv_NDdeconvolve-N( @g, @f ) );</lang>
 
Output:
2,392

edits