Re: Implementation of convn

[John W. Eaton]

On 23-Mar-2008, Søren Hauberg wrote:

| søn, 23 03 2008 kl. 22:35 +0100, skrev Søren Hauberg:
| > Hi,
| >   Attached is a fairly simple implementation of 'convn'. The
| > implementation consists of two files '__convn__.cc' and 'convn.m'. The
| > '__convn__' function computes the valid part of the convolution, while
| > the 'convn' function (the user interface) pads the data with zeros,
| > depending on the 'shape' argument. This split keeps the code fairly
| > simple, but it isn't as fast as humanly possible (fast enough though
| > IMHO). Is this approach acceptable?
| >   I should mention that the current code is missing
| >     * documentation
| >     * some optimisations to handle when B is larger than A
| > I'll implement this later this evening, and resend the code.
| Okay, attached is the updated code.
| 
| Søren

Will you please update according to the following comments and submit
this as an hg changeset?

| ## -*- texinfo -*-
| ## @deftypefn {Function File} @var{C} = convn(@var{A}, @var{B}, @var{shape})

Please don't use upper-case variables names in Octave code.

| function C = convn(A, B, shape = "full")
|   ## Check input
|   if (nargin < 2)
|     error("convn: not enough input arguments");
|   endif

To match the style of the rest of Octave, please use a space before
the open paren for function argument lists, like this:

  function c = convn (a, b, shape = "full")

    ## Check input
    if (nargin < 2)
      error ("convn: not enough input arguments");
    endif

| #include <octave/oct.h>

Files in the Octave sources should not use oct.h.  Instead, they
should include config.h (conditionally; see the other .cc files in the
Octave sources) and only the header files they need.

| inline
| int MAX(const int a, const int b)
| {
|   if (a > b)
|     return  a;
|   else
|     return b;
| }

The standard C++ header <algorithm> provides

  template <typename T> std::max (const T&, const T&)

so you should not need to define your own.

| template <class MatrixType, class SumType> 

Although we have used some StudlyCaps names in liboctave, I wish I
hadn't, so please don't add more.

You should probably avoid the use of MatrixType for the name of a
template parameter because that is already the name of a class in
Octave and could easily lead to confusion (it tripped me up for a few
minutes).

| octave_value
| convn(const MatrixType &A, const MatrixType &B)
| {
|   // Get sizes
|   const octave_idx_type ndims = A.ndims();
|   const octave_idx_type B_numel = B.numel();
|   const dim_vector A_size = A.dims();
|   const dim_vector B_size = B.dims();
|   
|   // Check input
|   if (ndims != B.ndims())
|     {
|       error("__convn__: first and second argument must have same 
dimensionality");
|       return octave_value();
|     }
| 
|   // Allocate output
|   dim_vector out_size(A_size);
|   for (octave_idx_type n = 0; n < ndims; n++)
|     out_size(n) = MAX(A_size(n) - B_size(n) + 1, 0);
|   MatrixType out = MatrixType(out_size);
|   const octave_idx_type out_numel = out.numel();
|   // Iterate over every element of 'out'.
|   dim_vector idx_dim(ndims);
|   Array<octave_idx_type> A_idx(idx_dim);
|   Array<octave_idx_type> B_idx(idx_dim);
|   Array<octave_idx_type> out_idx(idx_dim, 0);
|   for (octave_idx_type i = 0; i < out_numel; i++)
|     {

Please consider using some vertical whitesapce to make this easier to
read.

|       // For each neighbour
|       SumType sum = 0;
|       for (octave_idx_type n = 0; n < ndims; n++)
|         B_idx(n) = 0;
|       for (octave_idx_type j = 0; j < B_numel; j++)
|         {
|           for (octave_idx_type n = 0; n < ndims; n++)
|             A_idx(n) = out_idx(n) + (B_size(n)-1-B_idx(n));
|           sum += A(A_idx)*B(B_idx);
|           B.increment_index(B_idx, B_size);
|       }
|             
|       // Compute filter result
|       out(out_idx) = sum;
| 
|       // Prepare for next iteration
|       out.increment_index(out_idx, out_size);
|       
|       OCTAVE_QUIT;
|     }
|     
|   return octave_value(out);
| }

Is there a way to do this without increment_index?  That function will
make this code slow.

I think the prevailing style is to use OCTAVE_QUIT at the top of
loops.

Please use a space before the open paren for function argument lists.

| DEFUN_DLD(__convn__, args, , "\
| -*- texinfo -*-\n\
| @deftypefn {Loadable Function} __convn__(@var{A}, @var{B})\n\
| N-dimensional convolution. Only the valid part is computed. This is an 
internal\n\

To match the style used in the rest of Octave, please write

  DEFUN_DLD (__convn__, args, ,
    "-*- texinfo -*-\n\
  @deftypefn {Loadable Function} {} __convn__ (@var{a}, @var{b})\n\
  ...")

Note the addition of {} to the @deftypefn arguments (it is for the
return values).  Note also that doc strings for internal functions
should just say (for example)

    "-*- texinfo -*-\n\
  @deftypefn {Loadable Function} {} __convn__ (@var{a}, @var{b})\n\
  Undocumented internal function"

Currently we have a mixture of styles for these, so maybe they should
all be fixed to be consistent.

| function, and should not be called directly. Use @code{convn} instead.\n\
| @seealso{convn}\n\
| @end deftypefn\n\
| ")
| {
|   octave_value_list retval;
|   if (args.length() != 2)
|     {
|       print_usage ();
|       return retval;
|     }
|     
|   // Take action depending on input type
|   if (args(0).is_real_matrix() && args(1).is_real_matrix())
|     {
|       const NDArray A = args(0).array_value();
|       const NDArray B = args(1).array_value();
|       retval(0) = convn<NDArray, double>(A, B);
|     }
|   else if (args(0).is_complex_matrix() && args(1).is_complex_matrix())
|     {
|       const ComplexNDArray A = args(0).complex_matrix_value();
|       const ComplexNDArray B = args(1).complex_matrix_value();
|       retval(0) = convn<ComplexNDArray, Complex>(A, B);
|     }
|   else
|     {
|       error("__convn__: first and second input should be real, or complex 
arrays");
|       return retval;
|     }
|     
|   return retval;
| }

Please consider using the style

  if (args.length () == 2)
    {
      if (real matrix case)
        ...
      else if (complex matrix case)
        ...
      else
        error (...);
    }

so that there is only one return from the function.

Thanks,

jwe