Re: proposed implementation of interpn

[John W. Eaton]

On 12-Feb-2007, Alexander Barth wrote:

| here is a revised implementation of interpn.

OK, I included this function in Octave with a few more changes.


  bool isvector(const NDArray array)

I think you want "const NDArray& array" here.


  //  lookup a value in a sorted table (lookup.m)
  int lookup(double* x, const int& n, const double& y)

There is rarely any need to pass scalar values by const reference, the
first argument should be const, and the index values should have the
type octave_idx_type, so I changed this to

  octave_idx_type
  lookup (const double *x, octave_idx_type n, double y)

In addition, I changed the int index values used in the function to
octave_idx_type.

Also, for Octave code, please write the declaration this way, with the
return type on a separate line and the function name at the beginning
of a line.  That way, it is easier to find function definitions using
"grep ^FUNCTION_NAME ...".


  // n-dimensional linear interpolation
  void lin_interpn(int& n, int* size, int* scale, int& Ni, double extrapval, 
double** x,double* v, double** y, double* vi) 

I think all of these except vi should be const, the scalar values
should not be const reference, and the index values should be
octave_idx_type, so I changed this to

  void
  lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale,
               octave_idx_type Ni, double extrapval, const double **x,
               const double *v, const double **y, double *vi)

I also changed other index variables used inside the function to
octave_idx_type.


  DEFUN_DLD (interpn, args, ,
             "-*- texinfo -*-\n\
  @deftypefn {Function File} address@hidden interpn(@var{x1}, @var{x2},..., 
@var{xn}, @var{v}, @var{y1}, @var{y2},..., @var{yn}) \n\

In Texinfo, please use @dots{} instead of "...".


The printed documentation will look better if you write

  @seealso{interp1, interp2, ndgrid}\n\

instead of

  @seealso{interp1,interp2,ndgrid}\n\


Since this function only returns a single value, you can write

  octave_value retval;
  ...
  retval = Vi

instead of

  octave_value_list retval;
  ....
  retval(0) = Vi


  if (nargin % 2 == 0)

I think you need to also intercept calls with 0 or 1 arguments,
otherwise calls like interpn() will crash Octave.


  OCTAVE_LOCAL_BUFFER (NDArray, X, n);
  OCTAVE_LOCAL_BUFFER (NDArray, Y, n);
  // x and y are C arrays for X and Y
  // C arrays appear to be faster than NDArrays
  OCTAVE_LOCAL_BUFFER (double*, x, n);
  OCTAVE_LOCAL_BUFFER (double*, y, n);
  OCTAVE_LOCAL_BUFFER (int, scale, n);
  OCTAVE_LOCAL_BUFFER (int, size, n);

  NDArray V = args(n).array_value();
  NDArray Vi = NDArray(args(n+1).dims());

  if (error_state)
    {
      print_usage ();
      return retval;
    }

I think x, y, and V should be const and scale and size should be
octave_idx_type.  Also, at this point, we should be able to give a
more informative error message instead of just "incorrect usage", but
I don't know precisely what the message should be.  Can you please
write a better one?


  double* v = V.fortran_vec();

Since this function doesn't write to the elements of v, you can write

  const double* v = V.data ();


  for (int i = 0; i < n; i++)
    {
      X[i] = args(i).array_value();      
      Y[i] = args(n+i+1).array_value();

      if (error_state)
        {
          print_usage ();
          return retval;
        }

I think we should also be able to give a better error message here.


      if ((~isvector(X[i])) && X[i].numel() != size[i])

I think you mean to write ! instead of ~ here.


          x[i] = X[i].fortran_vec();

Since we don't write to elements of x, I think this can be

          x[i] = X[i].data ();


A version with my changes is checked in.

It would also be good to have a few tests for this function.

Thanks,

jwe