Deconvolution/2D+: Difference between revisions

Line 1,354:

Quite frankly I'm fairly astonished that it actually works...

(be warned this contains an exciting mix of 0- and 1- based indexes)

<lang Phix>-- ~~demo\rosetta\Deconvolution.exw~~

-- demo\rosetta\Deconvolution.exw

with javascript_semantics

function m_size(sequence m)

function m_size(sequence m)

--

--

-- returns the size of a matrix as a list of lengths

--

--

sequence res = {}

sequence res = {}

object me = m

object me = m

while sequence(me) do

while sequence(me) do

res &= length(me)

res &= length(me)

me = me[1]

me = me[1]

end while

end while

return res

return res

end function

end function

function product(sequence s)

function make_coordset(sequence size)

--

--

-- multiply all elements of s together

-- returns all points in the matrix, in zero-based indexes,

--

-- eg {{0,0,0}..{3,3,5}} for a 4x4x6 matrix [96 in total]

integer res = s[1]

--

for i=2 to length(s) do

sequence res = {}

res *= s[i]

integer count = product(size)

end for

for i=0 to count-1 do

return res

sequence coords = {}

end function

integer j = i

for s=length(size) to 1 by -1 do

function make_coordset(sequence size)

integer dimension = size[s]

--

coords &= mod(j,dimension)

-- returns all points in the matrix, in zero-based indexes,

j = floor(j/dimension)

-- eg {{0,0,0}..{3,3,5}} for a 4x4x6 matrix [96 in total]

end for

--

coords = reverse(coords)

sequence res = {}

res = append(res,coords)

integer count = product(size)

end for

for i=0 to count-1 do

return res

sequence coords = {}

end function

integer j = i

for s=length(size) to 1 by -1 do

function row(sequence g, f, gs, gc, fs, hs)

integer dimension = size[s]

--

coords &= mod(j,dimension)

--# Assembles a row, which is one of the simultaneous equations that needs

j = floor(j/dimension)

--# to be solved by reducing the whole set to reduced row echelon form. Note

end for

--# that each row describes the equation for a single cell of the 'g' function.

coords = reverse(coords)

--#

res = append(res,coords)

--# Arguments:

end for

--# g The "result" matrix of the convolution being undone.

return res

--# h The known "input" matrix of the convolution being undone.

end function

--# gs The size descriptor of 'g', passed as argument for efficiency.

--# gc The coordinate in 'g' that we are generating the equation for.

function row(sequence g, f, gs, gc, fs, hs)

--# fs The size descriptor of 'f', passed as argument for efficiency.

--

--# ~~Assembles~~ a ~~row,~~ ~~which~~ ~~is one~~ of the ~~simultaneous~~ ~~equations~~ ~~that~~ ~~needs~~

--# hs The size descriptor of 'h' (the unknown "input" matrix), passed

--# as argument for efficiency.

--# 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.

sequence row = {},

--#

coords = make_coordset(hs)

--# Arguments:

for i=1 to length(coords) do

--# g The "result" matrix of the convolution being undone.

sequence hc = coords[i]

--# h The known "input" matrix of the convolution being undone.

object fn = f

--# gs The size descriptor of 'g', passed as argument for efficiency.

for k=1 to length(gc) do

--# gc The coordinate in 'g' that we are generating the equation for.

integer d = gc[k]-hc[k]

--# fs The size descriptor of 'f', passed as argument for efficiency.

if d<0 or d>=fs[k] then

--# hs The size descriptor of 'h' (the unknown "input" matrix), passed

fn = 0

--# as argument for efficiency.

exit

--

end if

sequence row = {},

fn = fn[d+1]

coords = make_coordset(hs)

end for

for i=1 to length(coords) do

row = append(row,fn)

sequence hc = coords[i]

end for

object fn = f

object gn = g

for k=1 to length(gc) do

for i=1 to length(gc) do

integer d = gc[k]-hc[k]

gn = gn[gc[i]+1]

if d<0 or d>=fs[k] then

end for

fn = 0

row = append(row,gn)

exit

return row

end if

end function

fn = fn[d+1]

end for

function toRREF(sequence m)

row = append(row,fn)

--

end for

-- [renamed] copy of Reduced_row_echelon_form.htm#Phix

object gn = g

-- plus one small tweak, as noted below, exit->return,

for i=1 to length(gc) do

-- not that said seems to make any actual difference.

gn = gn[gc[i]+1]

--

end for

integer lead = 1,

row = append(row,gn)

rows = length(m),

return row

cols = length(m[1])

end function

m = deep_copy(m)

for r=1 to rows do

function toRREF(sequence m)

if lead>=cols then exit end if

--

integer i = r

-- [renamed] copy of Reduced_row_echelon_form.htm#Phix

while m[i][lead]=0 do

-- plus one small tweak, as noted below, exit->return,

i += 1

-- not that said seems to make any actual difference.

if i=rows then

--

i = r

integer lead = 1,

lead += 1

rows = length(m),

-- if lead=cols then exit end if

cols = length(m[1])

if lead=cols then return m end if

for r=1 to rows do

end if

if lead>=cols then exit end if

end while

integer i = r

sequence mi = deep_copy(m[r]),

while m[i][lead]=0 do

mr = sq_div(m[i],m[i][lead])

i += 1

m[i] = mi

if i=rows then

m[r] = mr

i = r

for j=1 to rows do

lead += 1

if j!=r then

-- if lead=cols then exit end if

m[j] = sq_sub(m[j],sq_mul(m[j][lead],m[r]))

if lead=cols then return m end if

end if

end if

end for

end while

lead += 1

-- nb m[i] is assigned before m[r], which matters when i=r:

end for

{m[r],m[i]} = {sq_div(m[i],m[i][lead]),m[r]}

return m

for j=1 to rows do

end function

if j!=r then

m[j] = sq_sub(m[j],sq_mul(m[j][lead],m[r]))

function lset(sequence h, sequence idx, object v)

end if

-- helper routine: store v somewhere deep inside h

end for

h = deep_copy(h)

lead += 1

integer i1 = idx[1]+1

end for

if length(idx)=1 then

return m

h[i1] = v

end function

else

h[i1] = lset(h[i1],idx[2..$],v)

function lset(sequence h, sequence idx, object v)

end if

-- helper routine: store v somewhere deep inside h

return h

integer i1 = idx[1]+1

end function

if length(idx)=1 then

h[i1] = v

function deconvolve(sequence g, f)

else

--

h[i1] = lset(h[i1],idx[2..$],v)

--# Deconvolve a pair of matrixes. Solves for 'h' such that 'g = f convolve h'.

end if

--#

return h

--# Arguments:

end function

--# g The matrix of data to be deconvolved.

--# f The matrix describing the convolution to be removed.

function deconvolve(sequence g, f)

--

-- Compute the sizes of the various matrixes involved.

--# Deconvolve a pair of matrixes. Solves for 'h' such that 'g = f convolve h'.

sequence gsize = m_size(g),

--#

fsize = m_size(f),

--# Arguments:

hsize = sq_add(sq_sub(gsize,fsize),1)

--# g The matrix of data to be deconvolved.

--# f The matrix describing the convolution to be removed.

-- Prepare the set of simultaneous equations to solve

--

sequence toSolve = {},

-- Compute the sizes of the various matrixes involved.

coords = make_coordset(gsize)

sequence gsize = m_size(g),

for i=1 to length(coords) do

fsize = m_size(f),

toSolve = append(toSolve,row(g,f,gsize,coords[i],fsize,hsize))

hsize = sq_add(sq_sub(gsize,fsize),1)

end for

-- Prepare the set of simultaneous equations to solve

-- Solve the equations

sequence toSolve = {},

sequence solved = toRREF(toSolve)

coords = make_coordset(gsize)

for i=1 to length(coords) do

-- Create a result matrix of the right size

toSolve = append(toSolve,row(g,f,gsize,coords[i],fsize,hsize))

object h = 0

end for

for i=length(hsize) to 1 by -1 do

h = repeat(h,hsize[i])

-- Solve the equations

end for

sequence solved = toRREF(toSolve)

-- Fill the results from the equations into the result matrix

-- Create a result matrix of the right size

coords = make_coordset(hsize)

object h = 0

for i=1 to length(coords) do

for i=length(hsize) to 1 by -1 do

h = lset(h,coords[i],solved[i][$])

h = repeat(h,hsize[i])

end for

end for

return h

end function

-- Fill the results from the equations into the result matrix

coords = make_coordset(hsize)

constant f1 = { 6, -9, -7, -5},

for i=1 to length(coords) do

g1 = {-48, 84, -16, 95, 125, -70, 7, 29, 54, 10},

h = lset(h,coords[i],solved[i][$])

h1 = {-8, 2, -9, -2, 9, -8, -2}

end for

return h

if deconvolve(g1, f1)!=h1 then ?9/0 end if

end function

if deconvolve(g1, h1)!=f1 then ?9/0 end if

constant f1 = { 6, -9, -7, -5},

constant f2 = {{-5, 2,-2,-6,-7},

g1 = {-48, 84, -16, 95, 125, -70, 7, 29, 54, 10},

{ 9, 7,-6, 5,-7},

h1 = {-8, 2, -9, -2, 9, -8, -2}

{ 1,-1, 9, 2,-7},

{ 5, 9,-9, 2,-5},

if deconvolve(g1, f1)!=h1 then ?9/0 end if

{-8, 5,-2, 8, 5}},

if deconvolve(g1, h1)!=f1 then ?9/0 end if

g2 = {{ 40, -21, 53, 42, 105, 1, 87, 60, 39, -28},

{-92, -64, 19,-167, -71, -47, 128,-109, 40, -21},

constant f2 = {{-5, 2,-2,-6,-7},

{ 58, 85,-93, 37, 101, -14, 5, 37, -76, -56},

{ 9, 7,-6, 5,-7},

{-90,-135, 60,-125, 68, 53, 223, 4, -36, -48},

{ 1,-1, 9, 2,-7},

{ 78, 16, 7,-199, 156,-162, 29, 28,-103, -10},

{ 5, 9,-9, 2,-5},

{-62, -89, 69, -61, 66, 193, -61, 71, -8, -30},

{-8, 5,-2, 8, 5}},

{ 48, -6, 21, -9,-150, -22, -56, 32, 85, 25}},

g2 = {{ 40, -21, 53, 42, 105, 1, 87, 60, 39, -28},

h2 = {{-8, 1,-7,-2,-9, 4},

{-92, -64, 19,-167, -71, -47, 128,-109, 40, -21},

{ 4, 5,-5, 2, 7,-1},

{ 58, 85,-93, 37, 101, -14, 5, 37, -76, -56},

{-6,-3,-3,-6, 9, 5}}

{-90,-135, 60,-125, 68, 53, 223, 4, -36, -48},

{ 78, 16, 7,-199, 156,-162, 29, 28,-103, -10},

if deconvolve(g2, f2)!=h2 then ?9/0 end if

{-62, -89, 69, -61, 66, 193, -61, 71, -8, -30},

if deconvolve(g2, h2)!=f2 then ?9/0 end if

{ 48, -6, 21, -9,-150, -22, -56, 32, 85, 25}},

h2 = {{-8, 1,-7,-2,-9, 4},

constant f3 = {{{-9, 5, -8}, { 3, 5, 1}},

{ 4, 5,-5, 2, 7,-1},

{{-1, -7, 2}, {-5, -6, 6}},

{-6,-3,-3,-6, 9, 5}}

{{ 8, 5, 8}, {-2, -6, -4}}},

g3 = {{{ 54, 42, 53, -42, 85, -72},

if deconvolve(g2, f2)!=h2 then ?9/0 end if

{ 45,-170, 94, -36, 48, 73},

if deconvolve(g2, h2)!=f2 then ?9/0 end if

{-39, 65,-112, -16, -78, -72},

{ 6, -11, -6, 62, 49, 8}},

constant f3 = {{{-9, 5, -8}, { 3, 5, 1}},

{{-57, 49, -23, 52,-135, 66},

{{-1, -7, 2}, {-5, -6, 6}},

{-23, 127, -58, -5,-118, 64},

{{ 8, 5, 8}, {-2, -6, -4}}},

{ 87, -16, 121, 23, -41, -12},

g3 = {{{ 54, 42, 53, -42, 85, -72},

{-19, 29, 35,-148, -11, 45}},

{ 45,-170, 94, -36, 48, 73},

{{-55,-147,-146, -31, 55, 60},

{-39, 65,-112, -16, -78, -72},

{-88, -45, -28, 46, -26,-144},

{ 6, -11, -6, 62, 49, 8}},

{-12,-107, -34, 150, 249, 66},

{{-57, 49, -23, 52,-135, 66},

{ 11, -15, -34, 27, -78, -50}},

{-23, 127, -58, -5,-118, 64},

{{ 56, 67, 108, 4, 2, -48},

{ 87, -16, 121, 23, -41, -12},

{ 58, 67, 89, 32, 32, -8},

{-19, 29, 35,-148, -11, 45}},

{-42, -31,-103, -30, -23, -8},

{{-55,-147,-146, -31, 55, 60},

{ 6, 4, -26, -10, 26, 12}}},

{-88, -45, -28, 46, -26,-144},

h3 = {{{ -6, -8, -5, 9},

{-12,-107, -34, 150, 249, 66},

{ -7, 9, -6, -8},

{ 11, -15, -34, 27, -78, -50}},

{ 2, -7, 9, 8}},

{{ 56, 67, 108, 4, 2, -48},

{{ 7, 4, 4, -6},

{ 58, 67, 89, 32, 32, -8},

{ 9, 9, 4, -4},

{-42, -31,-103, -30, -23, -8},

{ -3, 7, -2, -3}}}

{ 6, 4, -26, -10, 26, 12}}},

h3 = {{{ -6, -8, -5, 9},

if deconvolve(g3, f3)!=h3 then ?9/0 end if

{ -7, 9, -6, -8},

if deconvolve(g3, h3)!=f3 then ?9/0 end if

{ 2, -7, 9, 8}},

{{ 7, 4, 4, -6},

ppOpt({pp_Nest,2,pp_IntFmt,"%3d"})

{ 9, 9, 4, -4},

pp(deconvolve(g3, f3))

{ -3, 7, -2, -3}}}

pp(deconvolve(g3, h3))

if deconvolve(g3, f3)!=h3 then ?9/0 end if

if deconvolve(g3, h3)!=f3 then ?9/0 end if

ppOpt({pp_Nest,2,pp_IntFmt,"%3d"})

pp(deconvolve(g3, f3))

pp(deconvolve(g3, h3))</lang>

<pre>