Re: The conv2 mess

[Jordi Gutiérrez Hermoso]

2012/4/20 Jordi Gutiérrez Hermoso <address@hidden>:
> Can anyone understand wtf is going on in liboctave/oct-convn.cc ?

I think I can answer some of these now...

>  wtf is an inner or an outer convolution?

Inner convolution is the parts that get convolved when you use
"valid". Outer convolution is the rest when you call conv2 without a
shape parameter.

>  wtf are the variables passed to convolve_nd? What is "acd" and
>   "bcd"?

I still don't quite get this, but I think they're the combined
dimension vectors of the a, b, and c arrays.

>  htf is an n-dimensional convolution computed recursively? I
>   understand the 2d case is a base case, but how are the
>   higher-dimensional cases done?

I understand this now, but it's difficult to explain succinctly
without drawings and wooden blocks. ;-)

What I find odd is that the base case of his recursion is 2d
convolution instead of 1d, but I guess either makes sense.

>  wtf is part of this written as Fortran functions? Is 6 Fortran
>   copy-pasted files better than one C++ template?

The reason this is rewritten in Fortran seems to be simply in order to
be able to directly call BLAS gaxpy functions without going through
Octave's API. I'm not sure this is a worthwile reason, but it might be.

>  wtf was this rewritten completely, and if it was, why does it
>   still have Andy Adler's copyright statement, if it has none
>   of his code in place, except one cargo-culted Shape enum?

It actually does still have some of Andy's original structure,
vaguely. Some vestigial structures from Andy's original code remain
there by accident.

>  Even with the attached patch, wtf are A and B only approximately
>   equal?
>
>   x = rand(100);
>   y = ones(5);
>   A = conv2(x,y)(5:end-4,5:end-4);
>   B = conv2(x,y,"valid");

The issue is that the "valid" shape goes through the *conv2i functions
and without it, the computation goes through *conv2o functions. The
additions are performed in each case in a different order, and since
floating-point addition is non-associative, small differences creep
up.

I think this is surprising and undesirable behaviour. The "inner" and
"outer" convolutions should be performed exactly the same, except that
the "inner" convolution should simply stop computations earlier (i.e.
just not compute the edges that the "outer" convolution includes). I
think a proper fix of bug #34893 should address this issue completely.

> So htf do we fix this?

My proposal is this: get rid of the "inner" and "outer" separate code
paths in the *conv2.f source files. Annotate, if possible, the
convolve_nd function in liboctave/oct-convn.cc.

Can something be done about the six copy-pasted Fortran files? I hate
when I have to fix the same bug in 6 different source files. I would
prefer to rewrite this with C++ templates, but I don't know if the
apparent motive for the Fortran files (calling BLAS functions
directly) is reason enough to have copy-pasted code.

Thanks for listening to my ranting,
- Jordi G. H.