Deconvolution/2D+: Difference between revisions
Content added Content deleted
SqrtNegInf (talk | contribs) m (→{{header|Raku}}: sync rref/trim-system with other tasks, general tidying) |
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
||
Line 82: | Line 82: | ||
=={{header|C}}== |
=={{header|C}}== |
||
Very tedious code: unpacks 2D or 3D matrix into a vector with padding, do 1D FFT, then pack result back into matrix. |
Very tedious code: unpacks 2D or 3D matrix into a vector with padding, do 1D FFT, then pack result back into matrix. |
||
< |
<syntaxhighlight lang="c">#include <stdio.h> |
||
#include <stdlib.h> |
#include <stdlib.h> |
||
#include <math.h> |
#include <math.h> |
||
Line 299: | Line 299: | ||
printf("\n"); |
printf("\n"); |
||
} |
} |
||
}*/</ |
}*/</syntaxhighlight>Output<syntaxhighlight lang="text">deconv3(g, f): |
||
-6 -8 -5 9 |
-6 -8 -5 9 |
||
-7 9 -6 -8 |
-7 9 -6 -8 |
||
Line 316: | Line 316: | ||
8 5 8 |
8 5 8 |
||
-2 -6 -4 </ |
-2 -6 -4 </syntaxhighlight> |
||
=={{header|D}}== |
=={{header|D}}== |
||
< |
<syntaxhighlight lang="d">import std.stdio, std.conv, std.algorithm, std.numeric, std.range; |
||
class M(T) { |
class M(T) { |
||
Line 509: | Line 509: | ||
writeln(" f = ", f); |
writeln(" f = ", f); |
||
writeln("g deconv h = ", g.deconvolute(h)); |
writeln("g deconv h = ", g.deconvolute(h)); |
||
}</ |
}</syntaxhighlight> |
||
''todo(may be not :): pretty print & convert to normal D array'' |
''todo(may be not :): pretty print & convert to normal D array'' |
||
{{out}} |
{{out}} |
||
Line 520: | Line 520: | ||
=={{header|Go}}== |
=={{header|Go}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 711: | Line 711: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 740: | Line 740: | ||
Actually it is a matter of setting up the linear equations and then solving them. |
Actually it is a matter of setting up the linear equations and then solving them. |
||
'''Implementation''': < |
'''Implementation''': <syntaxhighlight lang="j">deconv3 =: 4 : 0 |
||
sz =. x >:@-&$ y NB. shape of z |
sz =. x >:@-&$ y NB. shape of z |
||
poi =. ,<"1 ($y) ,"0/&(,@i.) sz NB. pair of indexes |
poi =. ,<"1 ($y) ,"0/&(,@i.) sz NB. pair of indexes |
||
Line 748: | Line 748: | ||
sz $ 1e_8 round ({:"1 %. }:"1) T1 |
sz $ 1e_8 round ({:"1 %. }:"1) T1 |
||
) |
) |
||
round=: [ * <.@%~</ |
round=: [ * <.@%~</syntaxhighlight> |
||
'''Data''': < |
'''Data''': <syntaxhighlight lang="j">h1=: _8 2 _9 _2 9 _8 _2 |
||
f1=: 6 _9 _7 _5 |
f1=: 6 _9 _7 _5 |
||
g1=: _48 84 _16 95 125 _70 7 29 54 10 |
g1=: _48 84 _16 95 125 _70 7 29 54 10 |
||
Line 806: | Line 806: | ||
_42 _31 _103 _30 _23 _8 |
_42 _31 _103 _30 _23 _8 |
||
6 4 _26 _10 26 12 |
6 4 _26 _10 26 12 |
||
)</ |
)</syntaxhighlight> |
||
'''Tests''': < |
'''Tests''': <syntaxhighlight lang="j"> h1 -: g1 deconv3 f1 |
||
1 |
1 |
||
h2 -: g2 deconv3 f2 |
h2 -: g2 deconv3 f2 |
||
1 |
1 |
||
h3 -: g3 deconv3 f3 NB. -: checks for matching structure and data |
h3 -: g3 deconv3 f3 NB. -: checks for matching structure and data |
||
1</ |
1</syntaxhighlight> |
||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Julia has a deconv() function that works on Julia's builtin multidimensional arrays, but not on the nested type 2D and 3D arrays used in the task. So, the solution function, deconvn(), sets up repackaging for 1D fft. The actual solving work is done on one line of ifft/fft, and the rest of the code is merely to repackage the nested arrays. |
Julia has a deconv() function that works on Julia's builtin multidimensional arrays, but not on the nested type 2D and 3D arrays used in the task. So, the solution function, deconvn(), sets up repackaging for 1D fft. The actual solving work is done on one line of ifft/fft, and the rest of the code is merely to repackage the nested arrays. |
||
< |
<syntaxhighlight lang="julia">using FFTW, DSP |
||
const h1 = [-8, 2, -9, -2, 9, -8, -2] |
const h1 = [-8, 2, -9, -2, 9, -8, -2] |
||
Line 908: | Line 908: | ||
println(deconvn(f3nested, g3nested, |
println(deconvn(f3nested, g3nested, |
||
(length(h3nested), length(h3nested[1]), length(h3nested[1][1])))) # 3D |
(length(h3nested), length(h3nested[1]), length(h3nested[1][1])))) # 3D |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
Array{Array{Int64,1},1}[[[-6, -8, -5, 9], [-7, 9, -6, -8], [2, -7, 9, 8]], [[7, 4, 4, -6], [9, 9, 4, -4], [-3, 7, -2, -3]]] |
Array{Array{Int64,1},1}[[[-6, -8, -5, 9], [-7, 9, -6, -8], [2, -7, 9, 8]], [[7, 4, 4, -6], [9, 9, 4, -4], [-3, 7, -2, -3]]] |
||
Line 914: | Line 914: | ||
=={{header|Mathematica}} / {{header|Wolfram Language}}== |
=={{header|Mathematica}} / {{header|Wolfram Language}}== |
||
< |
<syntaxhighlight lang="mathematica">Round[ListDeconvolve[{6, -9, -7, -5}, {-48, 84, -16, 95, 125, -70, 7, 29, 54, 10}, Method -> "Wiener"]] |
||
Round[ListDeconvolve[{{-5, 2, -2, -6, -7}, {9, 7, -6, 5, -7}, {1, -1, 9, 2, -7}, {5, 9, -9, 2, -5}, {-8, 5, -2, 8, 5}}, |
Round[ListDeconvolve[{{-5, 2, -2, -6, -7}, {9, 7, -6, 5, -7}, {1, -1, 9, 2, -7}, {5, 9, -9, 2, -5}, {-8, 5, -2, 8, 5}}, |
||
Line 927: | Line 927: | ||
{-19, 29, 35, -148, -11, 45}}, {{-55, -147, -146, -31, 55, 60}, {-88, -45, -28, 46, -26, -144}, |
{-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}, |
{-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}}}, Method -> "Wiener"]]</ |
{-42, -31, -103, -30, -23, -8}, {6, 4, -26, -10, 26, 12}}}, Method -> "Wiener"]]</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 945: | Line 945: | ||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
{{trans|D}} |
{{trans|D}} |
||
< |
<syntaxhighlight lang="nim">import sequtils, typetraits |
||
type Size = uint64 |
type Size = uint64 |
||
Line 1,155: | Line 1,155: | ||
echo "f == g deconv h ? ", f == g.deconvolute(h) |
echo "f == g deconv h ? ", f == g.deconvolute(h) |
||
echo " f = ", f |
echo " f = ", f |
||
echo "g deconv f = ", g.deconvolute(h)</ |
echo "g deconv f = ", g.deconvolute(h)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,167: | Line 1,167: | ||
{{libheader|ntheory}} |
{{libheader|ntheory}} |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="perl">use feature 'say'; |
||
use ntheory qw/forsetproduct/; |
use ntheory qw/forsetproduct/; |
||
Line 1,317: | Line 1,317: | ||
pretty_print(0,@h); |
pretty_print(0,@h); |
||
print "\nff =\n"; |
print "\nff =\n"; |
||
pretty_print(0,@ff);</ |
pretty_print(0,@ff);</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>3D arrays: |
<pre>3D arrays: |
||
Line 1,354: | Line 1,354: | ||
Quite frankly I'm fairly astonished that it actually works...<br> |
Quite frankly I'm fairly astonished that it actually works...<br> |
||
(be warned this contains an exciting mix of 0- and 1- based indexes) |
(be warned this contains an exciting mix of 0- and 1- based indexes) |
||
<!--< |
<!--<syntaxhighlight lang="phix">(phixonline)--> |
||
<span style="color: #000080;font-style:italic;">-- demo\rosetta\Deconvolution.exw</span> |
<span style="color: #000080;font-style:italic;">-- demo\rosetta\Deconvolution.exw</span> |
||
<span style="color: #008080;">with</span> <span style="color: #000000;">javascript_semantics</span> |
<span style="color: #008080;">with</span> <span style="color: #000000;">javascript_semantics</span> |
||
Line 1,572: | Line 1,572: | ||
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f3</span><span style="color: #0000FF;">))</span> |
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f3</span><span style="color: #0000FF;">))</span> |
||
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">h3</span><span style="color: #0000FF;">))</span> |
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">h3</span><span style="color: #0000FF;">))</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,600: | Line 1,600: | ||
https://math.stackexchange.com/questions/380720/is-deconvolution-simply-division-in-frequency-domain |
https://math.stackexchange.com/questions/380720/is-deconvolution-simply-division-in-frequency-domain |
||
< |
<syntaxhighlight lang="python"> |
||
""" |
""" |
||
Line 1,731: | Line 1,731: | ||
pprint.pprint(deconv(g,f)) |
pprint.pprint(deconv(g,f)) |
||
</syntaxhighlight> |
|||
</lang> |
|||
Output: |
Output: |
||
Line 1,750: | Line 1,750: | ||
Translation of Tcl. |
Translation of Tcl. |
||
<lang |
<syntaxhighlight lang="raku" line># Deconvolution of N dimensional matrices. |
||
sub deconvolve-N ( @g, @f ) { |
sub deconvolve-N ( @g, @f ) { |
||
my @hsize = @g.shape »-« @f.shape »+» 1; |
my @hsize = @g.shape »-« @f.shape »+» 1; |
||
Line 1,888: | Line 1,888: | ||
my @ff = deconvolve-N( @g, @h-shaped ); |
my @ff = deconvolve-N( @g, @h-shaped ); |
||
say "\nff ="; |
say "\nff ="; |
||
pretty-print( @ff );</ |
pretty-print( @ff );</syntaxhighlight> |
||
Output: |
Output: |
||
Line 1,925: | Line 1,925: | ||
The trick to doing this (without using a library to do all the legwork for you) is to recast the higher-order solutions into solutions in the 1D case. This is done by regarding an ''n''-dimensional address as a coding of a 1-D address. |
The trick to doing this (without using a library to do all the legwork for you) is to recast the higher-order solutions into solutions in the 1D case. This is done by regarding an ''n''-dimensional address as a coding of a 1-D address. |
||
< |
<syntaxhighlight lang="tcl">package require Tcl 8.5 |
||
namespace path {::tcl::mathfunc ::tcl::mathop} |
namespace path {::tcl::mathfunc ::tcl::mathop} |
||
Line 2,089: | Line 2,089: | ||
return $h |
return $h |
||
}</ |
}</syntaxhighlight> |
||
Demonstrating how to use for the 3-D case: |
Demonstrating how to use for the 3-D case: |
||
< |
<syntaxhighlight lang="tcl"># A pretty-printer |
||
proc pretty matrix { |
proc pretty matrix { |
||
set size [rank $matrix] |
set size [rank $matrix] |
||
Line 2,142: | Line 2,142: | ||
# Now do the deconvolution and print it out |
# Now do the deconvolution and print it out |
||
puts h:\ [pretty [deconvolve $g $f]]</ |
puts h:\ [pretty [deconvolve $g $f]]</syntaxhighlight> |
||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 2,158: | Line 2,158: | ||
the [http://www.netlib.org/lapack/lug/node27.html <code>dgelsd</code>] function from the Lapack library. |
the [http://www.netlib.org/lapack/lug/node27.html <code>dgelsd</code>] function from the Lapack library. |
||
The <code>break</code> function breaks a long list into a sequence of sublists according to a given template, and the <code>band</code> function is taken from the [[Deconvolution/1D]] solution. |
The <code>break</code> function breaks a long list into a sequence of sublists according to a given template, and the <code>band</code> function is taken from the [[Deconvolution/1D]] solution. |
||
< |
<syntaxhighlight lang="ursala">#import std |
||
#import nat |
#import nat |
||
Line 2,171: | Line 2,171: | ||
lapack..dgelsd^^( |
lapack..dgelsd^^( |
||
(~&||0.!**+ ~&B^?a\~&Y@a ^lriFhNSS2iDrlYSK7LS2SL2rQ/~&alt band@alh2faltPrDPMX)^|\~&+ gang, |
(~&||0.!**+ ~&B^?a\~&Y@a ^lriFhNSS2iDrlYSK7LS2SL2rQ/~&alt band@alh2faltPrDPMX)^|\~&+ gang, |
||
@t =>~&l ~&L+@r))</ |
@t =>~&l ~&L+@r))</syntaxhighlight> |
||
The equations tend to become increasingly sparse in higher dimensions, |
The equations tend to become increasingly sparse in higher dimensions, |
||
so the following alternative implementation uses the sparse matrix |
so the following alternative implementation uses the sparse matrix |
||
Line 2,177: | Line 2,177: | ||
instead of Lapack, which is also callable in Ursala, adjusted as shown for the different |
instead of Lapack, which is also callable in Ursala, adjusted as shown for the different |
||
[http://www.basis.uklinux.net/avram/refman/umf-input-parameters.html calling convention]. |
[http://www.basis.uklinux.net/avram/refman/umf-input-parameters.html calling convention]. |
||
< |
<syntaxhighlight lang="ursala">deconv = # takes a number n to the n-dimensional deconvolution function |
||
~&?\math..div! iota; ~&!*; @h|\; -+ |
~&?\math..div! iota; ~&!*; @h|\; -+ |
||
Line 2,189: | Line 2,189: | ||
~&rFS+ num*rSS+ zipt^*D/~&r ^lrrSPT/~<K33tx zipt^/~&r ~&lSNyCK33+ zipp0, |
~&rFS+ num*rSS+ zipt^*D/~&r ^lrrSPT/~<K33tx zipt^/~&r ~&lSNyCK33+ zipp0, |
||
^/~&rx ~&B->NlNSPC ~&bt+-+-+-, |
^/~&rx ~&B->NlNSPC ~&bt+-+-+-, |
||
@t =>~&l ~&L+@r)+-</ |
@t =>~&l ~&L+@r)+-</syntaxhighlight> |
||
UMFPACK doesn't solve systems with more equations than unknowns, so |
UMFPACK doesn't solve systems with more equations than unknowns, so |
||
the system is pruned to a square matrix by first selecting an equation |
the system is pruned to a square matrix by first selecting an equation |
||
Line 2,200: | Line 2,200: | ||
However, some improvement may be possible by averaging the results over several runs. |
However, some improvement may be possible by averaging the results over several runs. |
||
Here is a test program. |
Here is a test program. |
||
< |
<syntaxhighlight lang="ursala">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.>>> |
||
f = <<<-9.,5.,-8.>,<3.,5.,1.>>,<<-1.,-7.,2.>,<-5.,-6.,6.>>,<<8.,5.,8.>,<-2.,-6.,-4.>>> |
f = <<<-9.,5.,-8.>,<3.,5.,1.>>,<<-1.,-7.,2.>,<-5.,-6.,6.>>,<<8.,5.,8.>,<-2.,-6.,-4.>>> |
||
Line 2,234: | Line 2,234: | ||
'h': deconv3(g,f), |
'h': deconv3(g,f), |
||
'f': deconv3(g,h)> |
'f': deconv3(g,h)> |
||
</syntaxhighlight> |
|||
</lang> |
|||
output: |
output: |
||
<pre style="height:25ex;overflow:scroll"> |
<pre style="height:25ex;overflow:scroll"> |
||
Line 2,287: | Line 2,287: | ||
{{libheader|Wren-complex}} |
{{libheader|Wren-complex}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
< |
<syntaxhighlight lang="ecmascript">import "/complex" for Complex |
||
import "/fmt" for Fmt |
import "/fmt" for Fmt |
||
Line 2,463: | Line 2,463: | ||
} |
} |
||
if (i < kx-1) System.print() |
if (i < kx-1) System.print() |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |