Deconvolution/2D+: Difference between revisions
m
→{{header|Wren}}: Minor tidy
m (→{{header|Wren}}: Minor tidy) |
|||
(16 intermediate revisions by 7 users not shown) | |||
Line 79:
[6, 4, -26, -10, 26, 12]]]
</pre>
=={{header|C}}==
Very tedious code: unpacks 2D or 3D matrix into a vector with padding, do 1D FFT, then pack result back into matrix.
<
#include <stdlib.h>
#include <math.h>
Line 300 ⟶ 299:
printf("\n");
}
}*/</
-6 -8 -5 9
-7 9 -6 -8
Line 317 ⟶ 316:
8 5 8
-2 -6 -4 </
=={{header|D}}==
<
class M(T) {
Line 510 ⟶ 509:
writeln(" f = ", f);
writeln("g deconv h = ", g.deconvolute(h));
}</
''todo(may be not :): pretty print & convert to normal D array''
{{out}}
Line 521 ⟶ 520:
=={{header|Go}}==
{{trans|C}}
<
import (
Line 712 ⟶ 711:
}
}
}</
{{out}}
Line 741 ⟶ 740:
Actually it is a matter of setting up the linear equations and then solving them.
'''Implementation''': <
sz =. x >:@-&$ y NB. shape of z
poi =. ,<"1 ($y) ,"0/&(,@i.) sz NB. pair of indexes
Line 749 ⟶ 748:
sz $ 1e_8 round ({:"1 %. }:"1) T1
)
round=: [ * <.@%~</
'''Data''': <
f1=: 6 _9 _7 _5
g1=: _48 84 _16 95 125 _70 7 29 54 10
Line 807 ⟶ 806:
_42 _31 _103 _30 _23 _8
6 4 _26 _10 26 12
)</
'''Tests''': <
1
h2 -: g2 deconv3 f2
1
h3 -: g3 deconv3 f3 NB. -: checks for matching structure and data
1</
=={{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
<
const h1 = [-8, 2, -9, -2, 9, -8, -2]
Line 868 ⟶ 866:
[-42, -31, -103, -30, -23, -8],
[6, 4, -26, -10, 26, 12]]]
function flatnested2d(a, siz)
ret = zeros(Int, prod(siz))
for i in 1:
ret[siz[2] * (i - 1) + j] = a[i][j]
end
Line 888 ⟶ 876:
function flatnested3d(a, siz)
ret = zeros(Int, prod(siz))
for i in 1:
ret[siz[2] * siz[3] * (i - 1) + siz[3] * (j - 1) + k] = a[i][j][k]
end
Line 898 ⟶ 884:
topow2(siz) = map(x -> nextpow(2, x), siz)
deconv1d(f1, g1) = Int.(round.(deconv(Float64.(g1), Float64.(f1))))
function deconv2d(f2, g2, xd2)
siz = topow2([length(g2), length(g2[1])])
h2 = Int.(round.(real.(ifft(fft(flatnested2d(g2, siz)) ./ fft(flatnested2d(f2, siz))))))
[[h2[siz[2] * (i - 1) + j] for j in 1:xd2[2]] for i in 1:xd2[1]]
end
function deconv3d(f3, g3, xd3)
siz = topow2([length(g3), length(g3[1]), length(g3[1][1])])
h3 = Int.(round.(real.(ifft(fft(flatnested3d(g3, siz)) ./ fft(flatnested3d(f3, siz))))))
[[[h3[siz[2] * siz[3] *(i - 1) + siz[3] * (j - 1) + k] for k in 1:xd3[3]]
for j in 1:xd3[2]] for i in 1:xd3[1]]
end
deconvn(f, g, tup=()) = length(tup) < 2 ? deconv1d(f, g) :
length(tup) == 2 ? deconv2d(f, g, tup) :
length(tup) == 3 ? deconv3d(f, g, tup) :
println("Array nesting > 3D not supported")
deconvn(f1, g1) # 1D
deconvn(f2nested, g2nested, (length(h2nested), length(h2nested[1]))) # 2D
println(
(length(h3nested), length(h3nested[1]), length(h3nested[1][1])))) # 3D
</
<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]]]
</pre>
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<
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 965 ⟶ 927:
{-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}}}, Method -> "Wiener"]]</
{{out}}
Line 980 ⟶ 942:
{2, -7, 9, 8, 0, 0}, {0, 0, 0, 0, 0, 0}}, {{7, 4, 4, -6, 0, 0}, {9, 9, 4, -4, 0, 0}, {-3, 7, -2, -3, 0, 0},
{0, 0, 0, 0, 0, 0}}}</pre>
=={{header|Nim}}==
{{trans|D}}
<syntaxhighlight lang="nim">import sequtils, typetraits
type Size = uint64
type M[T: SomeNumber] = object
dims: seq[Size]
subsizes: seq[Size]
data: seq[T]
####################################################################################################
# Miscellaneous.
func dotProduct[T: SomeNumber](a, b: openArray[T]): T =
assert a.len == b.len
for i in 0..a.high:
result += a[i] * b[i]
####################################################################################################
# Operations on M objects.
func setDimensions(m: var M; dimensions: varargs[Size]) =
for dim in dimensions:
if dim == 0:
raise newException(IndexDefect, "wrong dimension: 0")
m.dims = @dimensions
m.subsizes = m.dims
for i in 0..dimensions.high:
m.subsizes[i] = m.dims[(i+1)..^1].foldl(a * b, Size(1))
let dlength = m.dims[0] * m.subsizes[0]
if Size(m.data.len) != dlength:
m.data.setLen(dlength)
#---------------------------------------------------------------------------------------------------
func initM(m: var M; dimensions: varargs[Size]) =
m.setDimensions(dimensions)
#---------------------------------------------------------------------------------------------------
func set1DArray(m: var M; t: varargs[m.T]) =
let minLen = min(m.data.len, t.len)
m.data.setLen(minLen)
m.data[0..<minLen] = t[0..<minLen]
#---------------------------------------------------------------------------------------------------
func seqToIdx(m: M; s: Size): seq[Size] =
var acc = s
for subsize in m.subsizes:
result.add(acc div subsize)
acc = acc mod subsize
#---------------------------------------------------------------------------------------------------
template size(m: M): Size = Size(m.data.len)
#---------------------------------------------------------------------------------------------------
func checkBounds(m: M; indexes: varargs[Size]): bool =
if indexes.len > m.dims.len:
return false
for i, dim in indexes:
if dim >= m.dims[i]:
return false
result = true
#---------------------------------------------------------------------------------------------------
func `[]`(m: M; indexes: varargs[Size]): m.T =
if not m.checkBounds(indexes):
raise newException(IndexDefect, "index out of range: " & $indexes)
m.data[dotProduct(indexes, m.subsizes)]
#---------------------------------------------------------------------------------------------------
func `[]=`(m: M; indexes: varargs[int]; val: m.T) =
if not m.checkBounds(indexes):
raise newException(IndexDefect, "index out of range: " & $indexes)
m.data[dotProduct(indexes, m.subsizes)] = val
#---------------------------------------------------------------------------------------------------
func `==`(a, b: M): bool = a.dims == b.dims and a.data == b.data
#---------------------------------------------------------------------------------------------------
func `$`(m: M): string = $m.data
####################################################################################################
# Convolution/deconvolution.
func convolute(h, f: M): M =
## Result is "g".
var dims = h.dims
for i in 0..dims.high:
dims[i] += f.dims[i] - 1
result.initM(dims)
let bound = result.size
for i in 0..<h.size:
let hIndexes = h.seqToIdx(i)
for j in 0..<f.size:
let fIndexes = f.seqToIdx(j)
for k in 0..dims.high:
dims[k] = hIndexes[k] + fIndexes[k]
let idx1d = dotProduct(dims, result.subsizes)
if idx1d < bound:
result.data[idx1d] += h.data[i] * f.data[j]
else:
break # Bound reached.
#---------------------------------------------------------------------------------------------------
func deconvolute(g, f: M): M =
## Result is "h".
var dims = g.dims
for i, d in dims:
if d + 1 <= f.dims[i]:
raise newException(IndexDefect, "a dimension is zero or negative")
for i in 0..dims.high:
dims[i] -= f.dims[i] - 1
result.initM(dims)
for i in 0..<result.size:
let iIndexes = result.seqToIdx(i)
result.data[i] = g[iIndexes]
for j in 0..<i:
let jIndexes = result.seqToIdx(j)
for k in 0..dims.high:
dims[k] = iIndexes[k] - jIndexes[k]
if f.checkBounds(dims):
result.data[i] -= result.data[j] * f[dims]
when result.T is SomeInteger:
result.data[i] = result.data[i] div f.data[0]
else:
result.data[i] /= f.data[0]
####################################################################################################
# Transformation of a sequence into an M object.
func fold[T](a: seq[T]; d: var seq[Size]): auto =
if d.len == 0:
d.add(Size(a.len))
when a.elementType is seq:
if a.len == 0:
raise newException(ValueError, "empty dimension")
d.add(Size(a[0].len))
for elem in a:
if elem.len != a[0].len:
raise newException(ValueError, "not rectangular")
result = fold(a.foldl(a & b), d)
else:
if Size(a.len) != d.foldl(a * b):
raise newException(ValueError, "not same size")
result = a
#---------------------------------------------------------------------------------------------------
func arrtoM[T](a: T): auto =
var dims: seq[Size]
let d = fold(a, dims)
var res: M[d.elementType]
res.initM(dims)
res.set1DArray(d)
return res
#———————————————————————————————————————————————————————————————————————————————————————————————————
const 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] ] ]
const F = @[ @[ @[-9, 5, -8], @[ 3, 5, 1] ],
@[ @[-1, -7, 2], @[-5, -6, 6] ],
@[ @[ 8, 5, 8], @[-2, -6, -4] ] ]
let h = arrToM(H)
let f = arrToM(F)
let g = h.convolute(f)
echo "g == f convolute h ? ", g == f.convolute(h)
echo "h == g deconv f ? ", h == g.deconvolute(f)
echo "f == g deconv h ? ", f == g.deconvolute(h)
echo " f = ", f
echo "g deconv f = ", g.deconvolute(h)</syntaxhighlight>
{{out}}
<pre>g == f convolute h ? true
h == g deconv f ? true
f == g deconv h ? true
f = @[-9, 5, -8, 3, 5, 1, -1, -7, 2, -5, -6, 6, 8, 5, 8, -2, -6, -4]
g deconv f = @[-9, 5, -8, 3, 5, 1, -1, -7, 2, -5, -6, 6, 8, 5, 8, -2, -6, -4]</pre>
=={{header|Perl}}==
{{libheader|ntheory}}
{{trans|
<
use ntheory qw/forsetproduct/;
Line 1,134 ⟶ 1,317:
pretty_print(0,@h);
print "\nff =\n";
pretty_print(0,@ff);</
{{out}}
<pre>3D arrays:
Line 1,167 ⟶ 1,350:
]</pre>
=={{header|
{{trans|Tcl}}
Quite frankly I'm fairly astonished that it actually works...<br>
(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: #008080;">with</span> <span style="color: #000000;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">m_size</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- returns the size of a matrix as a list of lengths
--</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
<span style="color: #004080;">object</span> <span style="color: #000000;">me</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m</span>
<span style="color: #008080;">while</span> <span style="color: #004080;">sequence</span><span style="color: #0000FF;">(</span><span style="color: #000000;">me</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">&=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">me</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">me</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">me</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">make_coordset</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">size</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- returns all points in the matrix, in zero-based indexes,
-- eg {{0,0,0}..{3,3,5}} for a 4x4x6 matrix [96 in total]
--</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">count</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">product</span><span style="color: #0000FF;">(</span><span style="color: #000000;">size</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">count</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">coords</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">j</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">i</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">size</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">dimension</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">size</span><span style="color: #0000FF;">[</span><span style="color: #000000;">s</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">coords</span> <span style="color: #0000FF;">&=</span> <span style="color: #7060A8;">mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">j</span><span style="color: #0000FF;">,</span><span style="color: #000000;">dimension</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">j</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">j</span><span style="color: #0000FF;">/</span><span style="color: #000000;">dimension</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">coords</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">reverse</span><span style="color: #0000FF;">(</span><span style="color: #000000;">coords</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">coords</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">row</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">g</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">gs</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">gc</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fs</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">hs</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
--# Assembles a row, which is one of the simultaneous equations that needs
--# to be solved by reducing the whole set to reduced row echelon form. Note
--# that each row describes the equation for a single cell of the 'g' function.
--#
--# Arguments:
--# g The "result" matrix of the convolution being undone.
--# h The known "input" matrix of the convolution being undone.
--# gs The size descriptor of 'g', passed as argument for efficiency.
--# gc The coordinate in 'g' that we are generating the equation for.
--# fs The size descriptor of 'f', passed as argument for efficiency.
--# hs The size descriptor of 'h' (the unknown "input" matrix), passed
--# as argument for efficiency.
--</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">row</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{},</span>
<span style="color: #000000;">coords</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">make_coordset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">hs</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">coords</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">hc</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">coords</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #004080;">object</span> <span style="color: #000000;">fn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gc</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">gc</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]-</span><span style="color: #000000;">hc</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">d</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">fs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">fn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">exit</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">fn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">fn</span><span style="color: #0000FF;">[</span><span style="color: #000000;">d</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">row</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">row</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fn</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #004080;">object</span> <span style="color: #000000;">gn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">g</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gc</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">gn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">gn</span><span style="color: #0000FF;">[</span><span style="color: #000000;">gc</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">row</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">row</span><span style="color: #0000FF;">,</span><span style="color: #000000;">gn</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">row</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">toRREF</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- [renamed] copy of Reduced_row_echelon_form.htm#Phix
-- plus one small tweak, as noted below, exit->return,
-- not that said seems to make any actual difference.
--</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">lead</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">rows</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">cols</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span>
<span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">rows</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">lead</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">cols</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">i</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">r</span>
<span style="color: #008080;">while</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">lead</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">0</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">i</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">rows</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">i</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">r</span>
<span style="color: #000000;">lead</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #000080;font-style:italic;">-- if lead=cols then exit end if</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">lead</span><span style="color: #0000FF;">=</span><span style="color: #000000;">cols</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">m</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">mi</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">]),</span>
<span style="color: #000000;">mr</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">lead</span><span style="color: #0000FF;">])</span>
<span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mi</span>
<span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mr</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">rows</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">r</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">][</span><span style="color: #000000;">lead</span><span style="color: #0000FF;">],</span><span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">]))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">lead</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">m</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">lset</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">h</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">idx</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">object</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- helper routine: store v somewhere deep inside h</span>
<span style="color: #000000;">h</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">h</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">i1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">idx</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]+</span><span style="color: #000000;">1</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">idx</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">h</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span>
<span style="color: #008080;">else</span>
<span style="color: #000000;">h</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">h</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i1</span><span style="color: #0000FF;">],</span><span style="color: #000000;">idx</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">..$],</span><span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">h</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">g</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
--# Deconvolve a pair of matrixes. Solves for 'h' such that 'g = f convolve h'.
--#
--# Arguments:
--# g The matrix of data to be deconvolved.
--# f The matrix describing the convolution to be removed.
--
-- Compute the sizes of the various matrixes involved.</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">gsize</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m_size</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">fsize</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m_size</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">hsize</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_add</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gsize</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fsize</span><span style="color: #0000FF;">),</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- Prepare the set of simultaneous equations to solve</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">toSolve</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{},</span>
<span style="color: #000000;">coords</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">make_coordset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gsize</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">coords</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">toSolve</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">toSolve</span><span style="color: #0000FF;">,</span><span style="color: #000000;">row</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">gsize</span><span style="color: #0000FF;">,</span><span style="color: #000000;">coords</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">fsize</span><span style="color: #0000FF;">,</span><span style="color: #000000;">hsize</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000080;font-style:italic;">-- Solve the equations</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">solved</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">toRREF</span><span style="color: #0000FF;">(</span><span style="color: #000000;">toSolve</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- Create a result matrix of the right size</span>
<span style="color: #004080;">object</span> <span style="color: #000000;">h</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">hsize</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">h</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">h</span><span style="color: #0000FF;">,</span><span style="color: #000000;">hsize</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000080;font-style:italic;">-- Fill the results from the equations into the result matrix</span>
<span style="color: #000000;">coords</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">make_coordset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">hsize</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">coords</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">h</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">h</span><span style="color: #0000FF;">,</span><span style="color: #000000;">coords</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">solved</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][$])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">h</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">f1</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">5</span><span style="color: #0000FF;">},</span>
<span style="color: #000000;">g1</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">48</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">84</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">16</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">95</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">125</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">70</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">54</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">10</span><span style="color: #0000FF;">},</span>
<span style="color: #000000;">h1</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f1</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">h1</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">h1</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">f1</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">f2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{-</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">7</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">7</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">7</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">5</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">}},</span>
<span style="color: #000000;">g2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span> <span style="color: #000000;">40</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">21</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">53</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">42</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">105</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">87</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">60</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">39</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">28</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">92</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">64</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">167</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">71</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">47</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">128</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">109</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">40</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">21</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">58</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">85</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">93</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">37</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">101</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">14</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">37</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">76</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">56</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">90</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">135</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">60</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">125</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">68</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">53</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">223</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">36</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">48</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">78</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">16</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">199</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">156</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">162</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">28</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">103</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">10</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">62</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">89</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">69</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">61</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">66</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">193</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">61</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">71</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">30</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">48</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">150</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">22</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">56</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">32</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">85</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">25</span><span style="color: #0000FF;">}},</span>
<span style="color: #000000;">h2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">7</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">}}</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f2</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">h2</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">h2</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">f2</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">f3</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{{-</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">}},</span>
<span style="color: #0000FF;">{{-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">}},</span>
<span style="color: #0000FF;">{{</span> <span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8</span><span style="color: #0000FF;">},</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">4</span><span style="color: #0000FF;">}}},</span>
<span style="color: #000000;">g3</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{{</span> <span style="color: #000000;">54</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">42</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">53</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">42</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">85</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">72</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">45</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">170</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">94</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">36</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">48</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">73</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">39</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">65</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">112</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">16</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">78</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">72</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">11</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">62</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">49</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8</span><span style="color: #0000FF;">}},</span>
<span style="color: #0000FF;">{{-</span><span style="color: #000000;">57</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">49</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">23</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">52</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">135</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">66</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">23</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">127</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">58</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">118</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">64</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">87</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">16</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">121</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">41</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">12</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">19</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">35</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">148</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">11</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">45</span><span style="color: #0000FF;">}},</span>
<span style="color: #0000FF;">{{-</span><span style="color: #000000;">55</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">147</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">146</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">31</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">55</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">60</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">88</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">45</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">28</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">46</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">26</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">144</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">12</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">107</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">34</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">150</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">249</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">66</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">11</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">15</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">34</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">27</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">78</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">50</span><span style="color: #0000FF;">}},</span>
<span style="color: #0000FF;">{{</span> <span style="color: #000000;">56</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">67</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">108</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">48</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">58</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">67</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">89</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">32</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">32</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{-</span><span style="color: #000000;">42</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">31</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">103</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">30</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">23</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">26</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">26</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">12</span><span style="color: #0000FF;">}}},</span>
<span style="color: #000000;">h3</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{{</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8</span><span style="color: #0000FF;">}},</span>
<span style="color: #0000FF;">{{</span> <span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">9</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">4</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">}}}</span>
<span style="color: #008080;">if</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: #000000;">h3</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">if</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: #000000;">f3</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #7060A8;">ppOpt</span><span style="color: #0000FF;">({</span><span style="color: #004600;">pp_Nest</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_IntFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%3d"</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>
<!--</syntaxhighlight>-->
{{out}}
<pre>
{{{ -6, -8, -5, 9},
{ -7, 9, -6, -8},
{ 2, -7, 9, 8}},
{{ 7, 4, 4, -6},
{ 9, 9, 4, -4},
{ -3, 7, -2, -3}}}
{{{ -9, 5, -8},
{ 3, 5, 1}},
{{ -1, -7, 2},
{ -5, -6, 6}},
{{ 8, 5, 8},
{ -2, -6, -4}}}
</pre>
The version shipped in demo\rosetta contains the full 5 test sets: note that 5D takes a minute or two to complete.
=={{header|Python}}==
Tested on all 5 test cases.
Blows up with divide by zero error on 4d deconv(g,f) because
the fft(f) returns 0 for a sample. This shows the limits of
doing a deconvolution with fft.
https://math.stackexchange.com/questions/380720/is-deconvolution-simply-division-in-frequency-domain
<syntaxhighlight lang="python">
"""
https://rosettacode.org/wiki/Deconvolution/2D%2B
Working on 3 dimensional example using test data from the
RC task.
Python fft:
https://docs.scipy.org/doc/numpy/reference/routines.fft.html
"""
import numpy
import pprint
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]]]
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]]]
def trim_zero_empty(x):
"""
Takes a structure that represents an n dimensional example.
For a 2 dimensional example it will be a list of lists.
For a 3 dimensional one it will be a list of list of lists.
etc.
Actually these are multidimensional numpy arrays but I was thinking
in terms of lists.
Returns the same structure without trailing zeros in the inner
lists and leaves out inner lists with all zeros.
"""
if len(x) > 0:
if type(x[0]) != numpy.ndarray:
# x is 1d array
return list(numpy.trim_zeros(x))
else:
# x is a multidimentional array
new_x = []
for l in x:
tl = trim_zero_empty(l)
if len(tl) > 0:
new_x.append(tl)
return new_x
else:
# x is empty list
return x
def deconv(a, b):
"""
Returns function c such that b * c = a.
https://en.wikipedia.org/wiki/Deconvolution
"""
# Convert larger polynomial using fft
ffta = numpy.fft.fftn(a)
# Get it's shape so fftn will expand
# smaller polynomial to fit.
ashape = numpy.shape(a)
# Convert smaller polynomial with fft
# using the shape of the larger one
fftb = numpy.fft.fftn(b,ashape)
# Divide the two in frequency domain
fftquotient = ffta / fftb
# Convert back to polynomial coefficients using ifft
# Should give c but with some small extra components
c = numpy.fft.ifftn(fftquotient)
# Get rid of imaginary part and round up to 6 decimals
# to get rid of small real components
trimmedc = numpy.around(numpy.real(c),decimals=6)
# Trim zeros and eliminate
# empty rows of coefficients
cleanc = trim_zero_empty(trimmedc)
return cleanc
print("deconv(g,h)=")
pprint.pprint(deconv(g,h))
print(" ")
print("deconv(g,f)=")
pprint.pprint(deconv(g,f))
</syntaxhighlight>
Output:
<pre>
deconv(g,h)=
[[[-9.0, 5.0, -8.0], [3.0, 5.0, 1.0]],
[[-1.0, -7.0, 2.0], [-5.0, -6.0, 6.0]],
[[8.0, 5.0, 8.0], [-2.0, -6.0, -4.0]]]
deconv(g,f)=
[[[-6.0, -8.0, -5.0, 9.0], [-7.0, 9.0, -6.0, -8.0], [2.0, -7.0, 9.0, 8.0]],
[[7.0, 4.0, 4.0, -6.0], [9.0, 9.0, 4.0, -4.0], [-3.0, 7.0, -2.0, -3.0]]]
</pre>
=={{header|Raku}}==
(formerly Perl 6)
Translation of Tcl.
<syntaxhighlight lang="raku"
sub deconvolve-N ( @g, @f ) {
my @hsize = @g.shape »-« @f.shape »+» 1;
Line 1,184 ⟶ 1,763:
@h.AT-POS(|$_) = $v;
}
}
Line 1,195 ⟶ 1,774:
for ^@hc -> $i {
my $window = @gcoord[$i] - @hc[$i];
@fcoord.push($window) and next if 0
last;
}
Line 1,201 ⟶ 1,780:
}
@row.push: @g.AT-POS(|@gcoord);
}
# Constructs an AoA of coordinates to all elements of N dimensional array
sub coords ( @dim ) {
@[reverse $_ for [X] ([^$_] for reverse @dim)]
}
Line 1,212 ⟶ 1,791:
# Can handle over-specified systems (N unknowns in N + M equations)
sub rref (@m) {
my ($lead, $rows, $cols) = 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;
for ^$rows -> $n {
next if $n == $r;
@m[$n] »-=» @m[$r] »
}
++$lead;
}
}
return @m
my (\vars, @t) =
for
@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
if @array[0] ~~ Array {
say ' ' x $indent,"[";
say ' ' x $indent, "]{$indent??','!!''}";
} else {
Line 1,312 ⟶ 1,884:
my @h = deconvolve-N( @g, @f );
say "h =";
my @h-shaped[2;3;4] = @(deconvolve-N( @g, @f ));
my @ff = deconvolve-N( @g, @h-shaped );
say "\nff =";
Output:
Line 1,349 ⟶ 1,921:
],
]</pre>
=={{header|Tcl}}==
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.
<
namespace path {::tcl::mathfunc ::tcl::mathop}
Line 1,762 ⟶ 2,089:
return $h
}</
Demonstrating how to use for the 3-D case:
<
proc pretty matrix {
set size [rank $matrix]
Line 1,815 ⟶ 2,142:
# Now do the deconvolution and print it out
puts h:\ [pretty [deconvolve $g $f]]</
Output:
<pre>
Line 1,831 ⟶ 2,158:
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.
<
#import nat
Line 1,844 ⟶ 2,171:
lapack..dgelsd^^(
(~&||0.!**+ ~&B^?a\~&Y@a ^lriFhNSS2iDrlYSK7LS2SL2rQ/~&alt band@alh2faltPrDPMX)^|\~&+ gang,
@t =>~&l ~&L+@r))</
The equations tend to become increasingly sparse in higher dimensions,
so the following alternative implementation uses the sparse matrix
Line 1,850 ⟶ 2,177:
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].
<
~&?\math..div! iota; ~&!*; @h|\; -+
Line 1,862 ⟶ 2,189:
~&rFS+ num*rSS+ zipt^*D/~&r ^lrrSPT/~<K33tx zipt^/~&r ~&lSNyCK33+ zipp0,
^/~&rx ~&B->NlNSPC ~&bt+-+-+-,
@t =>~&l ~&L+@r)+-</
UMFPACK doesn't solve systems with more equations than unknowns, so
the system is pruned to a square matrix by first selecting an equation
Line 1,873 ⟶ 2,200:
However, some improvement may be possible by averaging the results over several runs.
Here is a test program.
<
f = <<<-9.,5.,-8.>,<3.,5.,1.>>,<<-1.,-7.,2.>,<-5.,-6.,6.>>,<<8.,5.,8.>,<-2.,-6.,-4.>>>
Line 1,907 ⟶ 2,234:
'h': deconv3(g,f),
'f': deconv3(g,h)>
</syntaxhighlight>
output:
<pre style="height:25ex;overflow:scroll">
Line 1,954 ⟶ 2,281:
<8.000000e+00,5.000000e+00,8.000000e+00>,
<-2.000000e+00,-6.000000e+00,-4.000000e+00>>>>
</pre>
=={{header|Wren}}==
{{trans|Go}}
{{libheader|Wren-complex}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./complex" for Complex
import "./fmt" for Fmt
var fft2 // recursive
fft2 = Fn.new { |buf, out, n, step, start|
if (step < n) {
fft2.call(out, buf, n, step*2, start)
fft2.call(out, buf, n, step*2, start + step)
var j = 0
while (j < n) {
var t = (Complex.imagMinusOne * Num.pi * j / n).exp * out[j+step+start]
buf[(j/2).floor + start] = out[j+start] + t
buf[((j+n)/2).floor + start] = out[j+start] - t
j = j + 2 * step
}
}
}
var fft = Fn.new { |buf, n|
var out = List.filled(n, null)
for (i in 0...n) out[i] = buf[i]
fft2.call(buf, out, n, 1, 0)
}
/* pad list length to power of two */
var padTwo = Fn.new { |g, le, nsl|
var n = 1
var ns = nsl[0]
if (ns != 0) {
n = ns
} else {
while (n < le) n = n * 2
}
var buf = List.filled(n, Complex.zero)
for (i in 0...le) buf[i] = Complex.new(g[i], 0)
nsl[0] = n
return buf
}
var deconv = Fn.new { |g, lg, f, lf, out, rowLe|
var ns = 0
var nsl = [ns]
var g2 = padTwo.call(g, lg, nsl)
var f2 = padTwo.call(f, lf, nsl)
ns = nsl[0]
fft.call(g2, ns)
fft.call(f2, ns)
var h = List.filled(ns, null)
for (i in 0...ns) h[i] = g2[i] / f2[i]
fft.call(h, ns)
for (i in 0...ns) {
if (h[i].real.abs < 1e-10) h[i] = Complex.zero
}
var i = 0
while (i > lf-lg-rowLe) {
out[-i] = (h[(i+ns)%ns]/32).real
i = i - 1
}
}
var unpack2 = Fn.new { |m, rows, le, toLe|
var buf = List.filled(rows*toLe, 0)
for (i in 0...rows) {
for (j in 0...le) buf[i*toLe+j] = m[i][j]
}
return buf
}
var pack2 = Fn.new { |buf, rows, fromLe, toLe, out|
for (i in 0...rows) {
for (j in 0...toLe) out[i][j] = buf[i*fromLe+j] / 4
}
}
var deconv2 = Fn.new { |g, rowG, colG, f, rowF, colF, out|
var g2 = unpack2.call(g, rowG, colG, colG)
var f2 = unpack2.call(f, rowF, colF, colG)
var ff = List.filled((rowG-rowF+1)*colG, 0)
deconv.call(g2, rowG*colG, f2, rowF*colG, ff, colG)
pack2.call(ff, rowG-rowF+1, colG, colG-colF+1, out)
}
var unpack3 = Fn.new { |m, x, y, z, toY, toZ|
var buf = List.filled(x*toY*toZ, 0)
for (i in 0...x) {
for (j in 0...y) {
for (k in 0...z) {
buf[(i*toY+j)*toZ+k] = m[i][j][k]
}
}
}
return buf
}
var pack3 = Fn.new { |buf, x, y, z, toY, toZ, out|
for (i in 0...x) {
for (j in 0...toY) {
for (k in 0...toZ) {
out[i][j][k] = buf[(i*y+j)*z+k] / 4
}
}
}
}
var deconv3 = Fn.new { |g, gx, gy, gz, f, fx, fy, fz, out|
var g2 = unpack3.call(g, gx, gy, gz, gy, gz)
var f2 = unpack3.call(f, fx, fy, fz, gy, gz)
var ff = List.filled((gx-fx+1)*gy*gz, 0)
deconv.call(g2, gx*gy*gz, f2, fx*gy*gz, ff, gy*gz)
pack3.call(ff, gx-fx+1, gy, gz, gy-fy+1, gz-fz+1, out)
}
var f = [
[[-9, 5, -8], [ 3, 5, 1]],
[[-1, -7, 2], [-5, -6, 6]],
[[ 8, 5, 8], [-2, -6, -4]]
]
var fx = f.count
var fy = f[0].count
var fz = f[0][0].count
var 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]]
]
var gx = g.count
var gy = g[0].count
var gz = g[0][0].count
var 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]]
]
var hx = gx - fx + 1
var hy = gy - fy + 1
var hz = gz - fz + 1
var h2 = List.filled(hx, null)
for (i in 0...hx) {
h2[i] = List.filled(hy, 0)
for (j in 0...hy) h2[i][j] = List.filled(hz, 0)
}
deconv3.call(g, gx, gy, gz, f, fx, fy, fz, h2)
System.print("deconv3(g, f):\n")
for (i in 0...hx) {
for (j in 0...hy) {
for (k in 0...hz) Fmt.write("$9.6h ", h2[i][j][k])
System.print()
}
if (i < hx-1) System.print()
}
var kx = gx - hx + 1
var ky = gy - hy + 1
var kz = gz - hz + 1
var f2 = List.filled(kx, null)
for (i in 0...kx) {
f2[i] = List.filled(ky, 0)
for (j in 0...ky) f2[i][j] = List.filled(kz, 0)
}
deconv3.call(g, gx, gy, gz, h, hx, hy, hz, f2)
System.print("\ndeconv3(g, h):\n")
for (i in 0...kx) {
for (j in 0...ky) {
for (k in 0...kz) Fmt.write("$9.6h ", f2[i][j][k])
System.print()
}
if (i < kx-1) System.print()
}</syntaxhighlight>
{{out}}
<pre>
deconv3(g, f):
-6 -8 -5 9
-7 9 -6 -8
2 -7 9 8
7 4 4 -6
9 9 4 -4
-3 7 -2 -3
deconv3(g, h):
-9 5 -8
3 5 1
-1 -7 2
-5 -6 6
8 5 8
-2 -6 -4
</pre>
|