Re: More octave-forge functions!!!

[David Bateman]

David Bateman wrote:

>With the bug fix to assignN I further tested gradient and found that I
>needed a fix for NDArrays as in the attached. Please replace the
>previous patch with the attached one..
>
>D.
>
>
>  
>
Here is a third attempt at gradient.m that also works with vectors... I
hope this is the last version!! Please use patch in previous version
with gradient.m replaced with this version.

As my first message also seems to have been blocked I include the text
of my original message below.

D.

David Bateman wrote:

>Here is several more of the octave-forge functions that are in the core
>of matlab converted for the core of octave. In addition this a a small
>typo in spline's help fixed... There were also a couple of supporting
>functions that I converted. The converted functions are
>
>  bicubic gradient inputname interp1 interp2 interpft mat2str ndgrid
>quadl pcg pchip
>  pchip_deriv pcr plot3 polyarea __plt3__
>
>Note that NDArrays and gradient don't work at the moment due to the bug
>I just reported on address@hidden in NDArray indexed assignment.
>
>Note ndgrid was from the maintainers mailing list recently and I'm just
>duplicating it in this patch..
>
>I also included pcr as well as pcg, even though only pcg is a core
>function, as pcr is the only iterative sparse solver we currently have
>that doesn't require a PD matrix.
>
>What I did was basically texinfo'fy all help strings, convert to octaves
>coding style, check matlab compatibility and treat inconsistencies where
>needed. Some of the inconsistencies I treated include
>
>* Add NDArray capabilities to a number of functions including gradient
>and interp1
>* Add the 'pp' argument to interp1
>* Convert __plt3__ to use the NaN values to flag missing data to gnuplot
>to largely speed up plot3
>
>There are a couple of little points in the converted functions, mainly
>with interp2, where I didn't implement 2-D spline interpolation and
>there is a comment about the bicubic function not allowing general meshs
>that I didn't address.
>
>This basically finishings all of the functions that I consider
>relatively simple to convert from octave-forge for use in octave. The
>remaining functions I consider to fall into two groups. The first of
>these groups of function in octave-forge which are core matlab functions is
>
>  brighten contourf drawnow fill fill3 gtext imread imwrite
>  legend meshc pareto patch pcolor peaks pie psi quiver
>  scatter stem surf surfc text view zoom
>
>Essentially these are all graphics functions and the graphics handles
>stuff will either largely change these functions, or simplify them and
>so I don't think these should be treated until the graphics handles
>stuff is addressed. Though they should go on someones todo list.
>
>The second group of functions in octave-forge that should be octave are
>
>  convhull convhulln csvread csvwrite del2 delaunay delaunay3
>  delaunayn dlmread dlmwrite edit ellipj ellipke expm1 fminbnd
>  funm fzero ginput griddata meshc ode23 ode45 optimset rat rats
>  sound soundsc textread tsearch voronoi voronoin waitbar xlsread
>  xmlread xmlwrite
>
>There are many reasons why not to port these yet, but basically I
>believe all need either additional work or discussion.
>
>For convhull, convgulln, delaunay, etc, which all depend on qhull there
>is the question of whether we want to add a dependency to octave on
>qhull. Otherwise they are essentially easy to port.
>
>For csvread, dlmread, csvwrite, dlmwrite it seems to me there must be a
>much simpler way of implementing these efficiently in a script file
>without the need for an oct-file.
>
>del2 should be relatively simple to convert, using the example of
>gradient.m, but it needs significant changes to allow NDArrays and be
>compatible to the current matlab version
>
>For edit I suspect some of the code here is not compatible with MinGW
>and so some changes are inevitably needed.
>
>For ellipj, expm1, do we want to add a dependency on GSL. If so then
>there are many many other functions that should be converted to use GSL
>as well. The list is really quite long, and should be discussed if the
>choice to make octave dependent on GSL is taken.
>
>For fminbnd there are significant changes to add the Options flags for
>full matlab compatibility.
>
>funm is fairly simple to convert, though function handles and inline
>functions should be added. However, it seems to me that thfm should be
>combined with funm to make the precision of calculation of trig
>functions better. However how to do this, particularly with function
>handles and inline functions doesn't seem obvious to me. Should we just
>take funm as is? Should in that case thfm also be added?
>
>I can't remember why I add fzero to my too hard list. Someone want to
>comment?
>
>ginput depends on an oct-file grab, that is dependent on Xlib. There
>should be a better way of doing this that is more portable if it is to
>be included in octave.
>
>There is a comment in ode23/ode45 of subtle differences in the interface
>between the octave and matlab versions. Do we want versions with such
>subtle differences? Would someone like to fix it to be completely
>compatible?
>
>Why port optimset if fminunc isn't ported.
>
>To port rat/rats need to convert to C++ include in pr-output.cc to allow
>"format rat" to work. This is just work, and should be fairly direct.
>
>sound/soundsc seem to be very system dependent, what is the best way of
>making these completely portable?
>
>Is the term code in waitbar portable? Should it be rewritten to use the
>GUI when written rather than like it currently is?
>
>xlsread is a very crippled implementation, that only takes tab separated
>text files.
>
>xmlread/xmlwrite need to be taken care of by someone else (cf mail
>archive for xmltree and xml4mat)
>
>I'm looking for both comments and help porting the above. So don't
>hesitate to respond..
>
>Regards
>D.
>
>  
>

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

## Copyright (C) 2000  Kai Habel
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
## 02110-1301, USA.

## -*- texinfo -*-
## @deftypefn {Function File} address@hidden = } gradient (@var{M})
## @deftypefnx {Function File} address@hidden, @var{y}, @dots{}] = } gradient 
(@var{M})
## @deftypefnx {Function File} address@hidden = } gradient (@var{M}, @var{s})
## @deftypefnx {Function File} address@hidden = } gradient (@var{M}, @var{dx}, 
@var{dy}, @dots{})
##
## Calculates the gradient. @address@hidden = gradient (@var{M})}
## calculates the one dimensional gradient if @var{M} is a vector. If
## @var{M} is a matrix the gradient is calculated for each row.
##
## @address@hidden, @var{y}] = gradient (@var{M})} calculates the one
## dimensional gradient for each direction if @var{M} if @var{M} is a
## matrix. Additional return arguments can be use for multi-dimensional
## matrices.
##
## Spacing values between two points can be provided by the
## @var{dx}, @var{dy} or @var{h} parameters. If @var{h} is supplied it
## is assumed to be the spacing in all directions. Otherwise, seperate
## values of the spacing can be supplied by the @var{dx}, etc variables.
## A scalar value specifies an equidistant spacing, while a vector value
## can be used to specify a variable spacing. The length must match
## their respective dimension of @var{M}.
## 
## At boundary points a linear extrapolation is applied. Interior points
## are calculated with the first approximation of the numerical gradient
##
## @example
## y'(i) = 1/(x(i+1)-x(i-1)) *(y(i-1)-y(i+1)).
## @end example
## 
## @end deftypefn

## Author:  Kai Habel <address@hidden>
## Modified: David Bateman <address@hidden> Added NDArray support

function [varargout] = gradient (M, varargin)
  
  if ((nargin < 1))
    print_usage ()
  endif

  transposed = false;
  if (isvector (M))
    ## make a column vector
    transposed = (size(M,2) == 1);
    M = M(:)';
  endif

  nd = ndims (M);
  sz = size (M);
  if (nargin > 2 && nargin != nd + 1)
    print_usage ()
  endif
    
  d = cell(1,nd);
  if (nargin == 1)
    for i=1:nd
      d{i} = ones(sz(i), 1);
    endfor
  elseif (nargin == 2)
    if (isscalar (varargin{1}))
      for i=1:nd
        d{i} = varargin{1} * ones(sz(i), 1);
      endfor
    else
      for i=1:nd
        d{i} = varargin{1};
      endfor
    endif
  else
    for i=1:nd
      if (isscalar (varargin{1}))
        d{i} = varargin{i} * ones(sz(i), 1);
      else
        d{i} = varargin{i};
      endif
    endfor

    ## Why the hell did matlab decide to swap these two values?
    tmp = d{1};
    d{1} = d{2};
    d{2} = tmp;
  endif

  for i = 1:max(2,min(nd,nargout))
    mr = sz(i);
    mc = prod([sz(1:i-1),sz(i+1:nd)]);
    Y = zeros (size(M), class(M));

    if (mr > 1)
      ## top and bottom boundary
      Y(1, :) = diff (M(1:2, :)) / d{i}(1);
      Y(mr, :) = diff (M(mr - 1:mr, :)) / d{i}(mr - 1);
    endif

    if (mr > 2)
      ## interior points
      Y(2:mr-1, :) = (M(3:mr, :) .- M(1:mr - 2, :)) ./ ...
          kron (d{i}(1:mr - 2) .+ d{i}(2:mr - 1), ones(1, mc));
    endif
    varargout{i} = ipermute (Y, [i:nd,1:i-1]);
    M = permute (M, [2:nd,1]);
  endfor

  ## Why the hell did matlab decide to swap these two values?
  tmp = varargout{1};
  varargout{1} = varargout{2};
  varargout{2} = tmp;

  if (transposed)
    varargout{1} = varargout{1}.';
  endif
endfunction