Mercurial > octave-nkf
changeset 13750:5f71ab377898 gui
Merge with default again
| author | Jordi Gutiérrez Hermoso <jordigh@octave.org> |
|---|---|
| date | Tue, 25 Oct 2011 11:28:41 -0700 |
| parents | 510237e67c2b (current diff) 62d1f56b0be7 (diff) |
| children | 43ffcaee3fea |
| files | NEWS configure.ac |
| diffstat | 24 files changed, 749 insertions(+), 467 deletions(-) [+] |
line wrap: on
line diff
--- a/.hgtags Sat Oct 22 16:22:04 2011 -0700 +++ b/.hgtags Tue Oct 25 11:28:41 2011 -0700 @@ -52,3 +52,4 @@ 695141f1c05cf1b240592bdd18e7a1503bb2a539 ss-3-3-55 901d466ee55ac902a875ec0ade6f1eccef0841dc release-3-4-1 3666e8e6f96e6899b8306d6ea9614aadf0500d67 release-3-4-2 +b0e70a71647b671ebcfa7a79af1ae6d3c0f52065 release-3-4-3
--- a/NEWS Sat Oct 22 16:22:04 2011 -0700 +++ b/NEWS Tue Oct 25 11:28:41 2011 -0700 @@ -51,12 +51,23 @@ in pattern and in string. Note that Matlab documents this behavior but the implementation does *not* always follow it. + ** Geometry functions derived from Qhull (convhull, delaunay, voronoi) + have been revamped. The options passed to the underlying qhull command + have been changed for better results or for Matlab compatibility. + + convhull : Default options are "Qt" for 2D, 3D, 4D inputs + Default options are "Qt Qx" for 5D and higher + delaunay : Default options are "Qt Qbb Qc Qz" for 2D and 3D inputs + Default options are "Qt Qbb Qc Qx" for 4D and higher + voronoi : No default arguments + ** Matlab-compatible preference functions: addpref getpref ispref rmpref setpref ** Other miscellaneous new functions: + nthargout iscolumn issrow zscore
--- a/doc/interpreter/func.txi Sat Oct 22 16:22:04 2011 -0700 +++ b/doc/interpreter/func.txi Tue Oct 25 11:28:41 2011 -0700 @@ -289,6 +289,12 @@ dimensions, and it is often desirable to give the individual return values distinct names. +It is possible to use the @code{nthargout} function to obtain only some +of the return values or several at once in a cell array. @ref{Cell Array +Objects} + +@DOCSTRING(nthargout) + In addition to setting @code{nargin} each time a function is called, Octave also automatically initializes @code{nargout} to the number of values that are expected to be returned. This allows you to write
--- a/scripts/general/module.mk Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/general/module.mk Tue Oct 25 11:28:41 2011 -0700 @@ -54,6 +54,7 @@ general/logspace.m \ general/nargchk.m \ general/nargoutchk.m \ + general/nthargout.m \ general/nextpow2.m \ general/num2str.m \ general/pol2cart.m \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/nthargout.m Tue Oct 25 11:28:41 2011 -0700 @@ -0,0 +1,113 @@ +## Copyright (C) 2011 Jordi Gutiérrez Hermoso +## +## 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 3 of the License, 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, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} nthargout (@var{n}, @var{func}, @dots{}) +## @deftypefnx {Function File} {} nthargout (@var{n}, @var{ntot}, @var{func}, @dots{}) +## Return the @var{n}th output argument of function given by the +## function handle or string @var{func}. Any arguments after @var{func} +## are passed to @var{func}. The total number of arguments to call +## @var{func} with can be passed in @var{ntot}; by default @var{ntot} +## is @var{n}. The input @var{n} can also be a vector of indices of the +## output, in which case the output will be a cell array of the +## requested output arguments. +## +## The intended use @code{nthargout} is to avoid intermediate variables. +## For example, when finding the indices of the maximum entry of a +## matrix, the following two compositions of nthargout +## +## @example +## @group +## @var{m} = magic (5); +## cell2mat (nthargout ([1, 2], @@ind2sub, size(@var{m}), +## nthargout (2, @@max, @var{m}(:)))) +## @result{} 5 3 +## @end group +## @end example +## +## @noindent +## are completely equivalent to the following lines: +## +## @example +## @group +## @var{m} = magic(5); +## [~, idx] = max (@var{M}(:)); +## [i, j] = ind2sub (size (@var{m}), idx); +## [i, j] +## @result{} 5 3 +## @end group +## @end example +## +## It can also be helpful to have all output arguments in a single cell +## in the following manner: +## +## @example +## @var{USV} = nthargout ([1:3], @@svd, hilb (5)); +## @end example +## +## @seealso{nargin, nargout, varargin, varargout, isargout} +## @end deftypefn + +## Author: Jordi Gutiérrez Hermoso + +function out = nthargout (n, varargin) + if (nargin < 2) + print_usage (); + endif + + if (isa (varargin{1}, "function_handle") || ischar (varargin{1})) + ntot = max (n(:)); + func = varargin{1}; + args = varargin(2:end); + elseif (isnumeric (varargin{1}) + && (isa (varargin{2}, "function_handle") || ischar (varargin{2}))) + ntot = varargin{1}; + func = varargin{2}; + args = varargin(3:end); + else + print_usage (); + endif + + if (any (n != fix (n)) || ntot != fix (ntot) || any (n <= 0) || ntot <= 0) + error ("nthargout: N and NTOT must consist of positive integers") + endif + + outargs = cell (1, ntot); + + try + [outargs{:}] = feval (func, args{:}); + if (numel (n) > 1) + out = outargs(n); + else + out = outargs{n}; + endif + catch + err = lasterr (); + if (strfind ("some elements undefined in return list", err)) + error ("nthargout: Too many output arguments: %d", ntot); + else + error (err); + endif + end_try_catch + +endfunction + +%!shared m +%! m = magic (5); +%!assert (nthargout ([1, 2], @ind2sub, size(m), nthargout (2, @max, m(:))), {5,3}) +%!assert (nthargout (3, @find, m(m>20)), [23, 24, 25, 21, 22]')
--- a/scripts/geometry/convhull.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/geometry/convhull.m Tue Oct 25 11:28:41 2011 -0700 @@ -18,35 +18,44 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {@var{H} =} convhull (@var{x}, @var{y}) -## @deftypefnx {Function File} {@var{H} =} convhull (@var{x}, @var{y}, @var{opt}) -## Return the index vector to the points of the enclosing convex hull. The -## data points are defined by the x and y vectors. +## @deftypefnx {Function File} {@var{H} =} convhull (@var{x}, @var{y}, @var{options}) +## Compute the convex hull of the set of points defined by the +## vectors @var{x} and @var{y}. The hull @var{H} is an index vector into +## the set of points and specifies which points form the enclosing hull. ## -## A third optional argument, which must be a string, contains extra options -## passed to the underlying qhull command. See the documentation for the -## Qhull library for details. +## An optional third argument, which must be a string or cell array of strings, +## contains options passed to the underlying qhull command. +## See the documentation for the Qhull library for details +## @url{http://www.qhull.org/html/qh-quick.htm#options}. +## The default option is @code{@{"Qt"@}}. ## -## @seealso{delaunay, convhulln} +## If @var{options} is not present or @code{[]} then the default arguments are +## used. Otherwise, @var{options} replaces the default argument list. +## To append user options to the defaults it is necessary to repeat the +## default arguments in @var{options}. Use a null string to pass no arguments. +## +## @seealso{convhulln, delaunay, voronoi} ## @end deftypefn ## Author: Kai Habel <kai.habel@gmx.de> -function H = convhull (x, y, opt) +function H = convhull (x, y, options) if (nargin != 2 && nargin != 3) print_usage (); endif - if (isvector (x) && isvector (y) && length (x) == length (y)) - if (nargin == 2) - i = convhulln ([x(:), y(:)]); - elseif (ischar (opt) || iscell (opt)) - i = convhulln ([x(:), y(:)], opt); - else - error ("convhull: third argument must be a string or cell array of strings"); - endif + if (! (isvector (x) && isvector (y) && length (x) == length (y)) + && ! size_equal (x, y)) + error ("convhull: X and Y must be the same size"); + elseif (nargin == 3 && ! (ischar (options) || iscellstr (options))) + error ("convhull: OPTIONS must be a string or cell array of strings"); + endif + + if (nargin == 2) + i = convhulln ([x(:), y(:)]); else - error ("convhull: first two input arguments must be vectors of same size"); + i = convhulln ([x(:), y(:)], options); endif n = rows (i); @@ -71,16 +80,21 @@ endfor H(n + 1) = H(1); + endfunction -%!testif HAVE_QHULL -%! x = -3:0.5:3; -%! y = abs (sin (x)); -%! assert (convhull (x, y, {"s","Qci","Tcv","Pp"}), [1;7;13;12;11;10;4;3;2;1]) %!demo %! x = -3:0.05:3; %! y = abs (sin (x)); %! k = convhull (x, y); -%! plot (x(k),y(k),'r-',x,y,'b+'); +%! plot (x(k),y(k),"r-;convex hull;", x,y,"b+;points;"); %! axis ([-3.05, 3.05, -0.05, 1.05]); + +%!testif HAVE_QHULL +%! x = -3:0.5:3; +%! y = abs (sin (x)); +%! assert (convhull (x, y), [1;7;13;12;11;10;4;3;2;1]) + +%% FIXME: Need input validation tests +
--- a/scripts/geometry/delaunay.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/geometry/delaunay.m Tue Oct 25 11:28:41 2011 -0700 @@ -17,74 +17,103 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{tri} =} delaunay (@var{x}, @var{y}) -## @deftypefnx {Function File} {@var{tri} =} delaunay (@var{x}, @var{y}, @var{opt}) -## The return matrix of size [n, 3] contains a set triangles which are -## described by the indices to the data point x and y vector. -## The triangulation satisfies the Delaunay circum-circle criterion. -## No other data point is in the circum-circle of the defining triangle. +## @deftypefn {Function File} {} delaunay (@var{x}, @var{y}) +## @deftypefnx {Function File} {@var{tri} =} delaunay (@var{x}, @var{y}) +## @deftypefnx {Function File} {@var{tri} =} delaunay (@var{x}, @var{y}, @var{options}) +## Compute the Delaunay triangulation for a 2-D set of points. +## The return value @var{tri} is a set of triangles which satisfies the +## Delaunay circum-circle criterion, i.e., only a single data point from +## [@var{x}, @var{y}] is within the circum-circle of the defining triangle. +## +## The set of triangles @var{tri} is a matrix of size [n, 3]. Each +## row defines a triangle and the three columns are the three vertices +## of the triangle. The value of @code{@var{tri}(i,j)} is an index into +## @var{x} and @var{y} for the location of the j-th vertex of the i-th +## triangle. ## -## A third optional argument, which must be a string, contains extra options -## passed to the underlying qhull command. See the documentation for the -## Qhull library for details. +## An optional third argument, which must be a string or cell array of strings, +## contains options passed to the underlying qhull command. +## See the documentation for the Qhull library for details +## @url{http://www.qhull.org/html/qh-quick.htm#options}. +## The default options are @code{@{"Qt", "Qbb", "Qc", "Qz"@}}. +## +## If @var{options} is not present or @code{[]} then the default arguments are +## used. Otherwise, @var{options} replaces the default argument list. +## To append user options to the defaults it is necessary to repeat the +## default arguments in @var{options}. Use a null string to pass no arguments. +## +## If no output argument is specified the resulting Delaunay triangulation +## is plotted along with the original set of points. ## ## @example ## @group ## x = rand (1, 10); -## y = rand (size (x)); +## y = rand (1, 10); ## T = delaunay (x, y); -## X = [x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1))]; -## Y = [y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1))]; +## VX = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ]; +## VY = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ]; ## axis ([0,1,0,1]); -## plot (X, Y, "b", x, y, "r*"); +## plot (VX, VY, "b", x, y, "r*"); ## @end group ## @end example -## @seealso{voronoi, delaunay3, delaunayn} +## @seealso{delaunay3, delaunayn, convhull, voronoi} ## @end deftypefn ## Author: Kai Habel <kai.habel@gmx.de> -function ret = delaunay (x, y, opt) +function tri = delaunay (x, y, options) if (nargin != 2 && nargin != 3) print_usage (); endif - if ((isvector (x) && isvector (y) && length (x) == length (y)) - || size_equal (x, y)) - if (nargin == 2) - tri = delaunayn ([x(:), y(:)]); - elseif (ischar (opt) || iscellstr (opt)) - tri = delaunayn ([x(:), y(:)], opt); - else - error ("delaunay: third argument must be a string"); - endif + if (! (isvector (x) && isvector (y) && length (x) == length (y)) + && ! size_equal (x, y)) + error ("delaunay: X and Y must be the same size"); + elseif (nargin == 3 && ! (ischar (options) || iscellstr (options))) + error ("delaunay: OPTIONS must be a string or cell array of strings"); + endif + + if (nargin == 2) + T = delaunayn ([x(:), y(:)]); else - error ("delaunay: first two input arguments must be matrices of same size"); + T = delaunayn ([x(:), y(:)], options); endif if (nargout == 0) x = x(:).'; y = y(:).'; - X = [x(tri(:,1)); x(tri(:,2)); x(tri(:,3)); x(tri(:,1))]; - Y = [y(tri(:,1)); y(tri(:,2)); y(tri(:,3)); y(tri(:,1))]; - plot(X, Y, 'b', x, y, 'r*'); + VX = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ]; + VY = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ]; + plot (VX, VY, "b", x, y, "r*"); else - ret = tri; + tri = T; endif + endfunction + +%!demo +%! old_state = rand ("state"); +%! restore_state = onCleanup (@() rand ("state", old_state)); +%! rand ("state", 1); +%! x = rand (1,10); +%! y = rand (1,10); +%! T = delaunay (x,y); +%! VX = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ]; +%! VY = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ]; +%! axis ([0,1,0,1]); +%! plot (VX,VY,"b", x,y,"r*"); + +%!testif HAVE_QHULL +%! x = [-1, 0, 1, 0]; +%! y = [0, 1, 0, -1]; +%! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,4;2,3,4]); + %!testif HAVE_QHULL %! x = [-1, 0, 1, 0, 0]; %! y = [0, 1, 0, -1, 0]; -%! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,5;1,4,5;2,3,5;3,4,5]) +%! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,5;1,4,5;2,3,5;3,4,5]); -%!demo -%! rand ('state', 1); -%! x = rand(1,10); -%! y = rand(size(x)); -%! T = delaunay(x,y); -%! X = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ]; -%! Y = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ]; -%! axis([0,1,0,1]); -%! plot(X,Y,'b',x,y,'r*'); +%% FIXME: Need input validation tests +
--- a/scripts/geometry/delaunay3.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/geometry/delaunay3.m Tue Oct 25 11:28:41 2011 -0700 @@ -17,42 +17,61 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{T} =} delaunay3 (@var{x}, @var{y}, @var{z}) -## @deftypefnx {Function File} {@var{T} =} delaunay3 (@var{x}, @var{y}, @var{z}, @var{opt}) -## A matrix of size [n, 4] is returned. Each row contains a -## set of tetrahedron which are -## described by the indices to the data point vectors (x,y,z). +## @deftypefn {Function File} {@var{tetr} =} delaunay3 (@var{x}, @var{y}, @var{z}) +## @deftypefnx {Function File} {@var{tetr} =} delaunay3 (@var{x}, @var{y}, @var{z}, @var{options}) +## Compute the Delaunay triangulation for a 3-D set of points. +## The return value @var{tetr} is a set of tetrahedrons which satisfies the +## Delaunay circum-circle criterion, i.e., only a single data point from +## [@var{x}, @var{y}, @var{z}] is within the circum-circle of the defining +## tetrahedron. ## -## A fourth optional argument, which must be a string or cell array of strings, -## contains extra options passed to the underlying qhull command. See the -## documentation for the Qhull library for details. -## @seealso{delaunay, delaunayn} +## The set of tetrahedrons @var{tetr} is a matrix of size [n, 4]. Each +## row defines a tetrahedron and the four columns are the four vertices +## of the tetrahedron. The value of @code{@var{tetr}(i,j)} is an index into +## @var{x}, @var{y}, @var{z} for the location of the j-th vertex of the i-th +## tetrahedron. +## +## An optional fourth argument, which must be a string or cell array of strings, +## contains options passed to the underlying qhull command. +## See the documentation for the Qhull library for details +## @url{http://www.qhull.org/html/qh-quick.htm#options}. +## The default options are @code{@{"Qt", "Qbb", "Qc", "Qz"@}}. +## +## If @var{options} is not present or @code{[]} then the default arguments are +## used. Otherwise, @var{options} replaces the default argument list. +## To append user options to the defaults it is necessary to repeat the +## default arguments in @var{options}. Use a null string to pass no arguments. +## +## @seealso{delaunay, delaunayn, convhull, voronoi} ## @end deftypefn ## Author: Kai Habel <kai.habel@gmx.de> -function tetr = delaunay3 (x, y, z, opt) +function tetr = delaunay3 (x, y, z, options) - if (nargin != 3 && nargin != 4) + if (nargin < 3 || nargin > 4) print_usage (); endif - if (isvector (x) && isvector (y) &&isvector (z) - && length (x) == length (y) && length(x) == length (z)) - if (nargin == 3) - tetr = delaunayn ([x(:), y(:), z(:)]); - elseif (ischar (opt) || iscell (opt)) - tetr = delaunayn ([x(:), y(:), z(:)], opt); - else - error ("delaunay3: fourth argument must be a string or cell array of strings"); - endif + if (! (isvector (x) && isvector (y) && isvector (z) + && length (x) == length (y) && length(x) == length (z))) + error ("delaunay: X, Y, and Z must be the same size"); + elseif (nargin == 4 && ! (ischar (options) || iscellstr (options))) + error ("delaunay3: OPTIONS must be a string or cell array of strings"); + endif + + if (nargin == 3) + tetr = delaunayn ([x(:), y(:), z(:)]); else - error ("delaunay3: first three input arguments must be vectors of same size"); + tetr = delaunayn ([x(:), y(:), z(:)], options); endif endfunction + %!testif HAVE_QHULL %! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1]; %! assert (sortrows (sort (delaunay3 (x, y, z), 2)), [1,2,3,4;1,2,4,5]) +%% FIXME: Need input validation tests +
--- a/scripts/geometry/delaunayn.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/geometry/delaunayn.m Tue Oct 25 11:28:41 2011 -0700 @@ -17,66 +17,77 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{t} =} delaunayn (@var{p}) -## @deftypefnx {Function File} {@var{t} =} delaunayn (@var{p}, @var{opt}) -## Form the Delaunay triangulation for a set of points. -## The Delaunay triangulation is a tessellation of the convex hull of the -## points such that no n-sphere defined by the n-triangles contains +## @deftypefn {Function File} {@var{T} =} delaunayn (@var{pts}) +## @deftypefnx {Function File} {@var{T} =} delaunayn (@var{pts}, @var{options}) +## Compute the Delaunay triangulation for an N-dimensional set of points. +## The Delaunay triangulation is a tessellation of the convex hull of a set +## of points such that no N-sphere defined by the N-triangles contains ## any other points from the set. -## The input matrix @var{p} of size @code{[n, dim]} contains @math{n} -## points in a space of dimension dim. The return matrix @var{t} has the -## size @code{[m, dim+1]}. It contains for each row a set of indices to -## the points, which describes a simplex of dimension dim. For example, -## a 2-D simplex is a triangle and 3-D simplex is a tetrahedron. +## +## The input matrix @var{pts} of size [n, dim] contains n points in a space of +## dimension dim. The return matrix @var{T} has size [m, dim+1]. Each row +## of @var{T} contains a set of indices back into the original set of points +## @var{pts} which describes a simplex of dimension dim. For example, a 2-D +## simplex is a triangle and 3-D simplex is a tetrahedron. ## -## Extra options for the underlying Qhull command can be specified by the -## second argument. This argument is a cell array of strings. The default -## options depend on the dimension of the input: +## An optional second argument, which must be a string or cell array of strings, +## contains options passed to the underlying qhull command. +## See the documentation for the Qhull library for details +## @url{http://www.qhull.org/html/qh-quick.htm#options}. +## The default options depend on the dimension of the input: ## ## @itemize -## @item 2D and 3D: @var{opt} = @code{@{"Qt", "Qbb", "Qc"@}} +## @item 2-D and 3-D: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}} ## -## @item 4D and higher: @var{opt} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}} +## @item 4-D and higher: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qx"@}} ## @end itemize ## -## If @var{opt} is [], then the default arguments are used. If @var{opt} -## is @code{@{"@w{}"@}}, then none of the default arguments are used by Qhull. -## See the Qhull documentation for the available options. +## If @var{options} is not present or @code{[]} then the default arguments are +## used. Otherwise, @var{options} replaces the default argument list. +## To append user options to the defaults it is necessary to repeat the +## default arguments in @var{options}. Use a null string to pass no arguments. ## -## All options can also be specified as single string, for example -## @code{"Qt Qbb Qc Qz"}. -## +## @seealso{delaunay, delaunay3, convhulln, voronoin} ## @end deftypefn -function t = delaunayn (p, varargin) +function T = delaunayn (pts, varargin) + if (nargin < 1) print_usage (); endif - t = __delaunayn__ (p, varargin{:}); + T = __delaunayn__ (pts, varargin{:}); - if (isa (p, "single")) + if (isa (pts, "single")) myeps = eps ("single"); else myeps = eps; endif - ## Try to remove the zero volume simplices. The volume of the i-th simplex is - ## given by abs(det(p(t(i,1:end-1),:)-p(t(i,2:end),:)))/prod(1:n) - ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a - ## relative volume less than some arbitrary criteria is rejected. The + ## Try to remove the zero volume simplices. The volume of the i-th simplex is + ## given by abs(det(pts(T(i,1:end-1),:)-pts(T(i,2:end),:)))/prod(1:n) + ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a + ## relative volume less than some arbitrary criteria is rejected. The ## criteria we use is the volume of the simplex corresponding to an ## orthogonal simplex is equal edge length all equal to the edge length of - ## the original simplex. If the relative volume is 1e3*eps then the simplex - ## is rejected. Note division of the two volumes means that the factor + ## the original simplex. If the relative volume is 1e3*eps then the simplex + ## is rejected. Note division of the two volumes means that the factor ## prod(1:n) is dropped. idx = []; - [nt, n] = size (t); + [nt, n] = size (T); + ## FIXME: Vectorize this for loop or convert to delaunayn to .oct function for i = 1:nt - X = p(t(i,1:end-1),:) - p(t(i,2:end),:); - if (abs (det (X)) / sqrt (sum (X .^ 2, 2)) < 1e3 * myeps) - idx = [idx, i]; + X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:); + if (abs (det (X)) / sqrt (sum (X .^ 2, 2)) < 1e3 * myeps) + idx(end+1) = i; endif endfor - t(idx,:) = []; + T(idx,:) = []; + endfunction + + +%% FIXME: Need tests for delaunayn + +%% FIXME: Need input validation tests +
--- a/scripts/geometry/griddata3.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/geometry/griddata3.m Tue Oct 25 11:28:41 2011 -0700 @@ -53,28 +53,32 @@ vi = griddatan ([x(:), y(:), z(:)], v(:), [xi(:), yi(:), zi(:)], varargin{:}); vi = reshape (vi, size (xi)); + endfunction + %!testif HAVE_QHULL -%! rand('state', 0); -%! x = 2 * rand(1000, 1) - 1; -%! y = 2 * rand(1000, 1) - 1; -%! z = 2 * rand(1000, 1) - 1; +%! old_state = rand ("state"); +%! restore_state = onCleanup (@() rand ("state", old_state)); +%! rand ("state", 0); +%! x = 2 * rand (1000, 1) - 1; +%! y = 2 * rand (1000, 1) - 1; +%! z = 2 * rand (1000, 1) - 1; %! v = x.^2 + y.^2 + z.^2; %! [xi, yi, zi] = meshgrid (-0.8:0.2:0.8); -%! ##vi = reshape (griddatan([x(:), y(:), z(:)], v, [xi(:), yi(:), zi(:)], 'linear'), size (xi)); %! vi = griddata3 (x, y, z, v, xi, yi, zi, 'linear'); %! vv = vi - xi.^2 - yi.^2 - zi.^2; -%! assert (max(abs(vv(:))), 0, 0.1) +%! assert (max (abs (vv(:))), 0, 0.1); %!testif HAVE_QHULL -%! rand('state', 0); -%! x = 2 * rand(1000, 1) - 1; -%! y = 2 * rand(1000, 1) - 1; -%! z = 2 * rand(1000, 1) - 1; +%! old_state = rand ("state"); +%! restore_state = onCleanup (@() rand ("state", old_state)); +%! rand ("state", 0); +%! x = 2 * rand (1000, 1) - 1; +%! y = 2 * rand (1000, 1) - 1; +%! z = 2 * rand (1000, 1) - 1; %! v = x.^2 + y.^2 + z.^2; %! [xi, yi, zi] = meshgrid (-0.8:0.2:0.8); -%! ##vi = reshape (griddatan([x(:), y(:), z(:)], v, [xi(:), yi(:), zi(:)], 'linear'), size (xi)); %! vi = griddata3 (x, y, z, v, xi, yi, zi, 'nearest'); %! vv = vi - xi.^2 - yi.^2 - zi.^2; -%! assert (max(abs(vv(:))), 0, 0.1) +%! assert (max (abs (vv(:))), 0, 0.1)
--- a/scripts/geometry/voronoi.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/geometry/voronoi.m Tue Oct 25 11:28:41 2011 -0700 @@ -18,18 +18,27 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} voronoi (@var{x}, @var{y}) -## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, "plotstyle") -## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, "plotstyle", @var{options}) +## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, @var{options}) +## @deftypefnx {Function File} {} voronoi (@dots{}, "linespec") +## @deftypefnx {Function File} {} voronoi (@var{hax}, @dots{}) +## @deftypefnx {Function File} {@var{h} =} voronoi (@dots{}) ## @deftypefnx {Function File} {[@var{vx}, @var{vy}] =} voronoi (@dots{}) -## Plot Voronoi diagram of points @code{(@var{x}, @var{y})}. +## Plot the Voronoi diagram of points @code{(@var{x}, @var{y})}. ## The Voronoi facets with points at infinity are not drawn. -## [@var{vx}, @var{vy}] = voronoi(@dots{}) returns the vertices instead of -## plotting the -## diagram. plot (@var{vx}, @var{vy}) shows the Voronoi diagram. +## +## If "linespec" is given it is used to set the color and line style of the +## plot. If an axis graphics handle @var{hax} is supplied then the Voronoi +## diagram is drawn on the specified axis rather than in a new figure. ## -## A fourth optional argument, which must be a string, contains extra options -## passed to the underlying qhull command. See the documentation for the -## Qhull library for details. +## The @var{options} argument, which must be a string or cell array of strings, +## contains options passed to the underlying qhull command. +## See the documentation for the Qhull library for details +## @url{http://www.qhull.org/html/qh-quick.htm#options}. +## +## If a single output argument is requested then the Voronoi diagram will be +## plotted and a graphics handle to the plot is returned. +## [@var{vx}, @var{vy}] = voronoi(@dots{}) returns the Voronoi vertices +## instead of plotting the diagram. ## ## @example ## @group @@ -57,23 +66,23 @@ ## Added optional fourth argument to pass options to the underlying ## qhull command -function [vvx, vvy] = voronoi (varargin) +function [vx, vy] = voronoi (varargin) if (nargin < 1) print_usage (); + elseif (nargout > 2) + error ("voronoi: No more than two output arguments supported"); endif narg = 1; if (isscalar (varargin{1}) && ishandle (varargin{1})) handl = varargin{1}; - narg++; if (! strcmp (get (handl, "type"), "axes")) error ("voronoi: expecting first argument to be an axes object"); endif - else - if (nargout < 2) - handl = gca (); - endif + narg++; + elseif (nargout < 2) + handl = gca (); endif if (nargin < 1 + narg || nargin > 3 + narg) @@ -87,24 +96,22 @@ if (narg <= nargin) if (iscell (varargin{narg})) opts = varargin(narg++); - elseif (ismatrix (varargin{narg})) - ## Accept but ignore the triangulation + elseif (isnumeric (varargin{narg})) + ## Accept, but ignore, the triangulation narg++; endif endif linespec = {"b"}; - if (narg <= nargin) - if (ischar (varargin{narg})) - linespec = varargin(narg); - endif + if (narg <= nargin && ischar (varargin{narg})) + linespec = varargin(narg); endif lx = length (x); ly = length (y); if (lx != ly) - error ("voronoi: arguments must be vectors of same length"); + error ("voronoi: X and Y must be vectors of the same length"); endif ## Add box to approximate rays to infinity. For Voronoi diagrams the @@ -126,12 +133,12 @@ [p, c, infi] = __voronoi__ ([[x(:) ; xbox(:)], [y(:); ybox(:)]], opts{:}); - idx = find (!infi); + idx = find (! infi); ll = length (idx); c = c(idx).'; k = sum (cellfun ("length", c)); - edges = cell2mat(cellfun (@(x) [x ; [x(end), x(1:end-1)]], c, - "uniformoutput", false)); + edges = cell2mat (cellfun (@(x) [x ; [x(end), x(1:end-1)]], c, + "uniformoutput", false)); ## Identify the unique edges of the Voronoi diagram edges = sortrows (sort (edges).').'; @@ -147,32 +154,34 @@ edges (:, edgeoutside) = []; ## Get points of the diagram - vx = reshape (p (edges, 1), size(edges)); - vy = reshape (p (edges, 2), size(edges)); + Vvx = reshape (p(edges, 1), size (edges)); + Vvy = reshape (p(edges, 2), size (edges)); if (nargout < 2) lim = [xmin, xmax, ymin, ymax]; - h = plot (handl, vx, vy, linespec{:}, x, y, '+'); + h = plot (handl, Vvx, Vvy, linespec{:}, x, y, '+'); axis (lim + 0.1 * [[-1, 1] * (lim (2) - lim (1)), ... [-1, 1] * (lim (4) - lim (3))]); if (nargout == 1) - vxx = h; + vx = h; endif - elseif (nargout == 2) - vvx = vx; - vvy = vy; else - error ("voronoi: only two or zero output arguments supported"); + vx = Vvx; + vy = Vvy; endif endfunction -%!testif HAVE_QHULL -%! phi=linspace(-pi,3/4*pi,8); -%! [x,y]=pol2cart(phi,1); -%! [vx,vy]=voronoi(x,y); -%! assert(vx(2,:),zeros(1,size(vx,2)),eps); -%! assert(vy(2,:),zeros(1,size(vy,2)),eps); %!demo %! voronoi (rand(10,1), rand(10,1)); + +%!testif HAVE_QHULL +%! phi = linspace (-pi, 3/4*pi, 8); +%! [x,y] = pol2cart (phi, 1); +%! [vx,vy] = voronoi (x,y); +%! assert(vx(2,:), zeros (1, columns (vx)), eps); +%! assert(vy(2,:), zeros (1, columns (vy)), eps); + +%% FIXME: Need input validation tests +
--- a/scripts/geometry/voronoin.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/geometry/voronoin.m Tue Oct 25 11:28:41 2011 -0700 @@ -20,14 +20,15 @@ ## @deftypefn {Function File} {[@var{C}, @var{F}] =} voronoin (@var{pts}) ## @deftypefnx {Function File} {[@var{C}, @var{F}] =} voronoin (@var{pts}, @var{options}) ## Compute N-dimensional Voronoi facets. The input matrix @var{pts} -## of size [n, dim] contains n points of dimension dim. +## of size [n, dim] contains n points in a space of dimension dim. ## @var{C} contains the points of the Voronoi facets. The list @var{F} -## contains for each facet the indices of the Voronoi points. +## contains, for each facet, the indices of the Voronoi points. ## -## A second optional argument, which must be a string, contains extra options -## passed to the underlying qhull command. See the documentation for the -## Qhull library for details. -## @seealso{voronoin, delaunay, convhull} +## An optional second argument, which must be a string or cell array of strings, +## contains options passed to the underlying qhull command. +## See the documentation for the Qhull library for details +## @url{http://www.qhull.org/html/qh-quick.htm#options}. +## @seealso{voronoi, convhulln, delaunayn} ## @end deftypefn ## Author: Kai Habel <kai.habel@gmx.de> @@ -43,17 +44,24 @@ print_usage (); endif - [np, dims] = size (pts); - if (np > dims) - if (nargin == 1) - [C, F, infi] = __voronoi__ (pts); - elseif (ischar (options) || iscellstr (options)) - [C, F, infi] = __voronoi__ (pts, options); - else - error ("voronoin: second argument must be a string or cell array of strings"); - endif + [np, dim] = size (pts); + + if (np <= dim) + error ("voronoin: number of points must be greater than their dimension"); + elseif (nargin == 2 && ! (ischar (options) || iscellstr (options))) + error ("voronoin: OPTIONS argument must be a string or cell array of strings"); + endif + if (nargin == 1) + [C, F] = __voronoi__ (pts); else - error ("voronoin: number of points must be greater than their dimension"); + [C, F] = __voronoi__ (pts, options); endif + endfunction + + +%% FIXME: Need functional tests + +%% FIXME: Need input validation tests +
--- a/scripts/linear-algebra/onenormest.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/linear-algebra/onenormest.m Tue Oct 25 11:28:41 2011 -0700 @@ -277,8 +277,10 @@ ## Only likely to be within a factor of 10. %!test +%! old_state = rand ("state"); +%! restore_state = onCleanup (@() rand ("state", old_state)); +%! rand ('state', 42); % Initialize to guarantee reproducible results %! N = 100; -%! rand ('state', 42); % Initialize to guarantee reproducible results %! A = rand (N); %! [nm1, v1, w1] = onenormest (A); %! [nminf, vinf, winf] = onenormest (A', 6);
--- a/scripts/plot/private/__axes_limits__.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/plot/private/__axes_limits__.m Tue Oct 25 11:28:41 2011 -0700 @@ -45,7 +45,7 @@ error ("%s: argument must be a 2 element vector", fcn); else if (arg(1) >= arg(2)) - error ("%s: axis limits must be increasing", fcn) + error ("%s: axis limits must be increasing", fcn); else set (h, fcn, arg(:)); endif
--- a/scripts/plot/trimesh.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/plot/trimesh.m Tue Oct 25 11:28:41 2011 -0700 @@ -39,28 +39,29 @@ triplot (tri, x, y, z, varargin{:}); else newplot (); - if (nargout > 0) - h = patch ("Vertices", [x(:), y(:), z(:)], "Faces", tri, - "FaceColor", "none", "EdgeColor", __next_line_color__(), - varargin{:}); - else - patch ("Vertices", [x(:), y(:), z(:)], "Faces", tri, - "FaceColor", "none", "EdgeColor", __next_line_color__(), - varargin{:}); - endif - + handle = patch ("Vertices", [x(:), y(:), z(:)], "Faces", tri, + "FaceColor", "none", "EdgeColor", __next_line_color__(), + varargin{:}); if (! ishold ()) set (gca(), "view", [-37.5, 30], "xgrid", "on", "ygrid", "on", "zgrid", "on"); endif + if (nargout > 0) + h = handle; + endif endif + endfunction + %!demo +%! old_state = rand ("state"); +%! restore_state = onCleanup (@() rand ("state", old_state)); +%! rand ("state", 10); %! N = 10; -%! rand ('state', 10) %! x = 3 - 6 * rand (N, N); %! y = 3 - 6 * rand (N, N); %! z = peaks (x, y); %! tri = delaunay (x(:), y(:)); %! trimesh (tri, x(:), y(:), z(:)); +
--- a/scripts/plot/triplot.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/plot/triplot.m Tue Oct 25 11:28:41 2011 -0700 @@ -22,7 +22,7 @@ ## @deftypefnx {Function File} {@var{h} =} triplot (@dots{}) ## Plot a triangular mesh in 2D@. The variable @var{tri} is the triangular ## meshing of the points @code{(@var{x}, @var{y})} which is returned from -## @code{delaunay}. If given, the @var{linespec} determines the properties +## @code{delaunay}. If given, @var{linespec} determines the properties ## to use for the lines. The output argument @var{h} is the graphic handle ## of the plot. ## @seealso{plot, trimesh, trisurf, delaunay} @@ -35,19 +35,24 @@ endif idx = tri(:, [1, 2, 3, 1]).'; - nt = size (tri, 1); + nt = rows (tri); + handle = plot ([x(idx); NaN(1, nt)](:), + [y(idx); NaN(1, nt)](:), varargin{:}); + if (nargout > 0) - h = plot ([x(idx); NaN(1, nt)](:), - [y(idx); NaN(1, nt)](:), varargin{:}); - else - plot ([x(idx); NaN(1, nt)](:), - [y(idx); NaN(1, nt)](:), varargin{:}); + h = handle; endif + endfunction + %!demo -%! rand ('state', 2) -%! x = rand (20, 1); -%! y = rand (20, 1); +%! old_state = rand ("state"); +%! restore_state = onCleanup (@() rand ("state", old_state)); +%! rand ("state", 2); +%! N = 20; +%! x = rand (N, 1); +%! y = rand (N, 1); %! tri = delaunay (x, y); %! triplot (tri, x, y); +
--- a/scripts/plot/trisurf.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/plot/trisurf.m Tue Oct 25 11:28:41 2011 -0700 @@ -27,7 +27,7 @@ ## @seealso{triplot, trimesh, delaunay3} ## @end deftypefn -function varargout = trisurf (tri, x, y, z, varargin) +function h = trisurf (tri, x, y, z, varargin) if (nargin < 3) print_usage (); @@ -55,11 +55,11 @@ varargin(end+(1:2)) = {"EdgeColor", "none"}; endif newplot (); - h = patch ("Faces", tri, "Vertices", [x(:), y(:), z(:)], - "FaceVertexCData", reshape (c, numel (c), 1), - varargin{:}); + handle = patch ("Faces", tri, "Vertices", [x(:), y(:), z(:)], + "FaceVertexCData", reshape (c, numel (c), 1), + varargin{:}); if (nargout > 0) - varargout = {h}; + h = handle; endif if (! ishold ()) @@ -67,11 +67,15 @@ "xgrid", "on", "ygrid", "on", "zgrid", "on"); endif endif + endfunction + %!demo +%! old_state = rand ("state"); +%! restore_state = onCleanup (@() rand ("state", old_state)); +%! rand ("state", 10); %! N = 10; -%! rand ('state', 10) %! x = 3 - 6 * rand (N, N); %! y = 3 - 6 * rand (N, N); %! z = peaks (x, y); @@ -99,4 +103,3 @@ %! tri = delaunay (x, y); %! trisurf (tri, x, y, z, "facecolor", "interp", "edgecolor", "k") -
--- a/scripts/sparse/svds.m Sat Oct 22 16:22:04 2011 -0700 +++ b/scripts/sparse/svds.m Tue Oct 25 11:28:41 2011 -0700 @@ -251,8 +251,12 @@ %! s = s(idx); %! u = u(:,idx); %! v = v(:,idx); +%! old_state1 = randn ("state"); +%! restore_state1 = onCleanup (@() randn ("state", old_state1)); +%! old_state2 = rand ("state"); +%! restore_state2 = onCleanup (@() rand ("state", old_state2)); %! randn ('state', 42); % Initialize to make normest function reproducible -%! rand ('state', 42) +%! rand ('state', 42); %! opts.v0 = rand (2*n,1); % Initialize eigs ARPACK starting vector %! % to guarantee reproducible results %!test @@ -279,3 +283,4 @@ %!test %! s = svds (speye (10)); %! assert (s, ones (6, 1), 2*eps); +
--- a/src/DLD-FUNCTIONS/__delaunayn__.cc Sat Oct 22 16:22:04 2011 -0700 +++ b/src/DLD-FUNCTIONS/__delaunayn__.cc Tue Oct 25 11:28:41 2011 -0700 @@ -62,8 +62,8 @@ DEFUN_DLD (__delaunayn__, args, , "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{T} =} __delaunayn__ (@var{P})\n\ -@deftypefnx {Loadable Function} {@var{T} =} __delaunayn__ (@var{P}, @var{opt})\n\ +@deftypefn {Loadable Function} {@var{T} =} __delaunayn__ (@var{pts})\n\ +@deftypefnx {Loadable Function} {@var{T} =} __delaunayn__ (@var{pts}, @var{options})\n\ Undocumented internal function.\n\ @end deftypefn") @@ -73,7 +73,6 @@ #ifdef HAVE_QHULL retval(0) = 0.0; - std::string options = ""; int nargin = args.length (); if (nargin < 1 || nargin > 2) @@ -86,67 +85,52 @@ const octave_idx_type dim = p.columns (); const octave_idx_type n = p.rows (); - // default options + // Default options + std::string options; if (dim <= 3) - options = "Qt Qbb Qc"; + options = "Qt Qbb Qc Qz"; else options = "Qt Qbb Qc Qx"; - if (nargin == 2) { - if (args(1).is_empty ()) - { - // keep default options - } - else if (args(1).is_string ()) - { - // option string is directly provided + if (args(1).is_string ()) options = args(1).string_value (); - } - else if (args(1).is_cell ()) - { - options = ""; - - Cell c = args(1).cell_value (); - for (octave_idx_type i = 0; i < c.numel (); i++) - { + else if (args(1).is_empty ()) + ; // Use default options + else if (args(1).is_cellstr ()) + { + options = ""; + Array<std::string> tmp = args(1).cellstr_value (); - if (! c.elem(i).is_string ()) - { - error ("__delaunayn__: all options must be strings"); - return retval; - } - - options = options + c.elem(i).string_value () + " "; - } - } - else - { - error ("__delaunayn__: OPT argument must be a string, cell array of strings, or empty"); - return retval; - } + for (octave_idx_type i = 0; i < tmp.numel (); i++) + options += tmp(i) + " "; + } + else + { + error ("__delaunayn__: OPTIONS argument must be a string, cell array of strings, or empty"); + return retval; + } } - //octave_stdout << "options " << options << std::endl; - if (n > dim + 1) { p = p.transpose (); double *pt_array = p.fortran_vec (); boolT ismalloc = false; - OCTAVE_LOCAL_BUFFER (char, flags, 250); + // Qhull flags argument is not const char* + OCTAVE_LOCAL_BUFFER (char, flags, 9 + options.length()); sprintf (flags, "qhull d %s", options.c_str ()); - // If you want some debugging information replace the 0 pointer - // with stdout or some other file open for writing. - + // Replace the 0 pointer with stdout for debugging information FILE *outfile = 0; FILE *errfile = stderr; - - if (! qh_new_qhull (dim, n, pt_array, ismalloc, flags, outfile, errfile)) + + int exitcode = qh_new_qhull (dim, n, pt_array, + ismalloc, flags, outfile, errfile); + if (! exitcode) { // triangulate non-simplicial facets qh_triangulate (); @@ -160,54 +144,48 @@ if (! facet->upperdelaunay) nf++; - // Double check + // Double check. Non-simplicial facets will cause segfault below if (! facet->simplicial) { error ("__delaunayn__: Qhull returned non-simplicial facets -- try delaunayn with different options"); + exitcode = 1; break; } } - Matrix simpl (nf, dim+1); - - FORALLfacets + if (! exitcode) { - if (! facet->upperdelaunay) - { - octave_idx_type j = 0; + Matrix simpl (nf, dim+1); - FOREACHvertex_ (facet->vertices) + FORALLfacets + { + if (! facet->upperdelaunay) { - // if delaunayn crashes, enable this check -#if 0 - if (j > dim) - { - error ("__delaunayn__: internal error. Qhull returned non-simplicial facets"); - return retval; - } -#endif + octave_idx_type j = 0; - simpl(i, j++) = 1 + qh_pointid(vertex->point); + FOREACHvertex_ (facet->vertices) + { + simpl(i, j++) = 1 + qh_pointid(vertex->point); + } + i++; } - i++; } - } - - retval(0) = simpl; - // free long memory - qh_freeqhull (! qh_ALL); - - // free short memory and memory allocator - int curlong, totlong; - qh_memfreeshort (&curlong, &totlong); - - if (curlong || totlong) - warning ("__delaunay__: did not free %d bytes of long memory (%d pieces)", - totlong, curlong); + retval(0) = simpl; + } } else error ("__delaunayn__: qhull failed"); + + // Free memory from Qhull + qh_freeqhull (! qh_ALL); + + int curlong, totlong; + qh_memfreeshort (&curlong, &totlong); + + if (curlong || totlong) + warning ("__delaunay__: did not free %d bytes of long memory (%d pieces)", + totlong, curlong); } else if (n == dim + 1) {
--- a/src/DLD-FUNCTIONS/__voronoi__.cc Sat Oct 22 16:22:04 2011 -0700 +++ b/src/DLD-FUNCTIONS/__voronoi__.cc Tue Oct 25 11:28:41 2011 -0700 @@ -53,10 +53,11 @@ #endif #endif -DEFUN_DLD (__voronoi__, args, , +DEFUN_DLD (__voronoi__, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{tri} =} __voronoi__ (@var{point})\n\ -@deftypefnx {Loadable Function} {@var{tri} =} __voronoi__ (@var{point}, @var{options})\n\ +@deftypefn {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{pts})\n\ +@deftypefnx {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{pts}, @var{options})\n\ +@deftypefnx {Loadable Function} {@var{C}, @var{F}, @var{Inf_Pts} =} __voronoi__ (@dots{})\n\ Undocumented internal function.\n\ @end deftypefn") { @@ -73,57 +74,50 @@ return retval; } - std::string options = "qhull v Fv T0"; + std::string options = ""; if (nargin == 2) { - if (args(1).is_cellstr ()) + if (args(1).is_string ()) + options = args(1).string_value (); + else if (args(1).is_empty ()) + ; // Use default options + else if (args(1).is_cellstr ()) { + options = ""; Array<std::string> tmp = args(1).cellstr_value (); for (octave_idx_type i = 0; i < tmp.numel (); i++) - options += " " + tmp(i); + options += tmp(i) + " "; } - else if (args(1).is_string ()) - options += " " + args(1).string_value (); else { - error ("__voronoi__: OPTIONS argument must be a string or cellstr"); + error ("__voronoi__: OPTIONS argument must be a string, cell array of strings, or empty"); return retval; } } Matrix p (args(0).matrix_value ()); - const octave_idx_type dim = p.columns (); const octave_idx_type np = p.rows (); + p = p.transpose (); - double *pt_array = p.fortran_vec (); - //double pt_array[dim * np]; - //for (int i = 0; i < np; i++) - // { - // for (int j = 0; j < dim; j++) - // { - // pt_array[j+i*dim] = p(i,j); - // } - // } - boolT ismalloc = false; - // If you want some debugging information replace the 0 pointer - // with stdout or some other file open for writing. - + // Replace the 0 pointer with stdout for debugging information FILE *outfile = 0; FILE *errfile = stderr; - // Qhull flags argument is not const char*... - OCTAVE_LOCAL_BUFFER (char, flags, options.length () + 1); + // Qhull flags argument is not const char* + OCTAVE_LOCAL_BUFFER (char, flags, 12 + options.length ()); - strcpy (flags, options.c_str ()); + sprintf (flags, "qhull v FV %s", options.c_str ()); - if (! qh_new_qhull (dim, np, pt_array, ismalloc, flags, outfile, errfile)) + int exitcode = qh_new_qhull (dim, np, pt_array, + ismalloc, flags, outfile, errfile); + if (! exitcode) { facetT *facet; vertexT *vertex; @@ -134,6 +128,7 @@ for (i = 0; i < np; i++) ni[i] = 0; qh_setvoronoi_all (); + bool infinity_seen = false; facetT *neighbor, **neighborp; coordT *voronoi_vertex; @@ -147,6 +142,7 @@ { if (qh hull_dim == 3) qh_order_vertexneighbors (vertex); + infinity_seen = false; FOREACHneighbor_ (vertex) @@ -173,9 +169,7 @@ for (octave_idx_type d = 0; d < dim; d++) v(0,d) = octave_Inf; - boolMatrix AtInf (np, 1); - for (i = 0; i < np; i++) - AtInf(i) = false; + boolMatrix AtInf (np, 1, false); octave_value_list F (np, octave_value ()); k = 0; i = 0; @@ -189,6 +183,7 @@ { if (qh hull_dim == 3) qh_order_vertexneighbors(vertex); + infinity_seen = false; RowVector facet_list (ni[k++]); m = 0; @@ -229,24 +224,27 @@ for (i = 0; i < r; i++) C.elem (i) = F(i); - retval(0) = v; + if (nargout == 3) + { + AtInf.resize (r, 1); + retval(2) = AtInf; + } retval(1) = C; - AtInf.resize (r, 1); - retval(2) = AtInf; - - // free long memory - qh_freeqhull (! qh_ALL); - - // free short memory and memory allocator - int curlong, totlong; - qh_memfreeshort (&curlong, &totlong); - - if (curlong || totlong) - warning ("__voronoi__: did not free %d bytes of long memory (%d pieces)", totlong, curlong); + retval(0) = v; } else error ("__voronoi__: qhull failed"); + // free memory from Qhull + qh_freeqhull (! qh_ALL); + + int curlong, totlong; + qh_memfreeshort (&curlong, &totlong); + + if (curlong || totlong) + warning ("__voronoi__: did not free %d bytes of long memory (%d pieces)", + totlong, curlong); + #else error ("__voronoi__: not available in this version of Octave"); #endif
--- a/src/DLD-FUNCTIONS/convhulln.cc Sat Oct 22 16:22:04 2011 -0700 +++ b/src/DLD-FUNCTIONS/convhulln.cc Tue Oct 25 11:28:41 2011 -0700 @@ -53,24 +53,39 @@ DEFUN_DLD (convhulln, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{h} =} convhulln (@var{p})\n\ -@deftypefnx {Loadable Function} {@var{h} =} convhulln (@var{p}, @var{opt})\n\ +@deftypefn {Loadable Function} {@var{h} =} convhulln (@var{pts})\n\ +@deftypefnx {Loadable Function} {@var{h} =} convhulln (@var{pts}, @var{options})\n\ @deftypefnx {Loadable Function} {[@var{h}, @var{v}] =} convhulln (@dots{})\n\ -Return an index vector to the points of the enclosing convex hull.\n\ -The input matrix of size [n, dim] contains n points of dimension dim.\n\n\ -If a second optional argument is given, it must be a string or cell array\n\ -of strings containing options for the underlying qhull command. (See\n\ -the Qhull documentation for the available options.) The default options\n\ -are \"s Qci Tcv\".\n\ -If the second output @var{v} is requested the volume of the convex hull is\n\ -calculated.\n\n\ -@seealso{convhull, delaunayn}\n\ +Compute the convex hull of the set of points @var{pts} which is a matrix\n\ +of size [n, dim] containing n points in a space of dimension dim.\n\ +The hull @var{h} is an index vector into the set of points and specifies\n\ +which points form the enclosing hull.\n\ +\n\ +An optional second argument, which must be a string or cell array of strings,\n\ +contains options passed to the underlying qhull command.\n\ +See the documentation for the Qhull library for details\n\ +@url{http://www.qhull.org/html/qh-quick.htm#options}.\n\ +The default options depend on the dimension of the input:\n\ +\n\ +@itemize\n\ +@item 2D, 3D, 4D: @var{options} = @code{@{\"Qt\"@}}\n\ +\n\ +@item 5D and higher: @var{options} = @code{@{\"Qt\", \"Qx\"@}}\n\ +@end itemize\n\ +\n\ +If @var{options} is not present or @code{[]} then the default arguments are\n\ +used. Otherwise, @var{options} replaces the default argument list.\n\ +To append user options to the defaults it is necessary to repeat the\n\ +default arguments in @var{options}. Use a null string to pass no arguments.\n\ +\n\ +If the second output @var{v} is requested the volume of the enclosing\n\ +convex hull is calculated.\n\n\ +@seealso{convhull, delaunayn, voronoin}\n\ @end deftypefn") { octave_value_list retval; #ifdef HAVE_QHULL - std::string options; int nargin = args.length (); if (nargin < 1 || nargin > 2) @@ -79,85 +94,88 @@ return retval; } + Matrix p (args(0).matrix_value ()); + const octave_idx_type dim = p.columns (); + const octave_idx_type n = p.rows (); + + // Default options + std::string options; + if (dim <= 4) + options = "Qt"; + else + options = "Qt Qx"; + if (nargin == 2) { - if (args (1).is_string ()) + if (args(1).is_string ()) options = args(1).string_value (); - else if (args(1).is_cell ()) + else if (args(1).is_empty ()) + ; // Use default options + else if (args(1).is_cellstr ()) { - Cell c = args(1).cell_value (); options = ""; - for (octave_idx_type i = 0; i < c.numel (); i++) - { - if (! c.elem(i).is_string ()) - { - error ("convhulln: OPT must be a string or cell array of strings"); - return retval; - } + Array<std::string> tmp = args(1).cellstr_value (); - options = options + c.elem(i).string_value() + " "; - } + for (octave_idx_type i = 0; i < tmp.numel (); i++) + options += tmp(i) + " "; } else { - error ("convhulln: OPT must be a string or cell array of strings"); + error ("convhulln: OPTIONS must be a string, cell array of strings, or empty"); return retval; } - } - else - // turn on some consistency checks - options = "s Qci Tcv"; + } - Matrix p (args(0).matrix_value ()); - const octave_idx_type dim = p.columns (); - const octave_idx_type n = p.rows (); p = p.transpose (); - double *pt_array = p.fortran_vec (); - - boolT ismalloc = False; - - std::ostringstream buf; + boolT ismalloc = false; - buf << "qhull QJ " << options; - - std::string buf_string = buf.str (); + // FIXME: we can't just pass options.c_str () to qh_new_qhull + // because the argument is not declared const. Ugh. Unless qh_new_qhull + // really needs to modify this argument, someone should fix QHULL. + OCTAVE_LOCAL_BUFFER (char, flags, 7 + options.length ()); - // FIXME -- we can't just pass buf_string.c_str () to qh_new_qhull - // because the argument is not declared const. Ugh. Unless - // qh_new_qhull really needs to modify this argument, someone should - // fix QHULL. + sprintf (flags, "qhull %s", options.c_str ()); - OCTAVE_LOCAL_BUFFER (char, flags, buf_string.length () + 1); - - strcpy (flags, buf_string.c_str ()); - - if (! qh_new_qhull (dim, n, pt_array, ismalloc, flags, 0, stderr)) + // Replace the 0 pointer with stdout for debugging information + FILE *outfile = 0; + FILE *errfile = stderr; + + int exitcode = qh_new_qhull (dim, n, pt_array, + ismalloc, flags, outfile, errfile); + if (! exitcode) { - // If you want some debugging information replace the NULL - // pointer with stdout - vertexT *vertex, **vertexp; facetT *facet; setT *vertices; - unsigned int nf = qh num_facets; + bool nonsimp_seen = false; + octave_idx_type nf = qh num_facets; - Matrix idx (nf, dim); + Matrix idx (nf, dim + 1); - octave_idx_type j, i = 0; + octave_idx_type i = 0, j; FORALLfacets { j = 0; - if (! facet->simplicial) - // should never happen with QJ - error ("convhulln: non-simplicial facet"); + + if (! nonsimp_seen && ! facet->simplicial) + { + nonsimp_seen = true; + if (options.find ("QJ") != std::string::npos) + { + // should never happen with QJ + error ("convhulln: qhull failed. Option 'QJ' returned non-simplicial facet"); + break; + } + } if (dim == 3) { vertices = qh_facet3vertex (facet); FOREACHvertex_ (vertices) idx(i, j++) = 1 + qh_pointid(vertex->point); + qh_settempfree (&vertices); } else @@ -174,11 +192,15 @@ } } if (j < dim) - // likewise but less fatal warning ("facet %d only has %d vertices", i, j); + i++; } + // Remove extra dimension if all facets were simplicial + if (! nonsimp_seen) + idx.resize (nf, dim, 0.0); + if (nargout == 2) // calculate volume of convex hull // taken from qhull src/geom2.c @@ -213,19 +235,21 @@ retval(1) = octave_value (qh totvol); } - retval(0) = octave_value (idx); + retval(0) = idx; } + else + error ("convhulln: qhull failed"); - // free long memory + // free memory from Qhull qh_freeqhull (! qh_ALL); - // free short memory and memory allocator int curlong, totlong; qh_memfreeshort (&curlong, &totlong); if (curlong || totlong) warning ("convhulln: did not free %d bytes of long memory (%d pieces)", - totlong, curlong); + totlong, curlong); + #else error ("convhulln: not available in this version of Octave"); #endif @@ -236,10 +260,23 @@ /* %!testif HAVE_QHULL %! cube = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1]; -%! [h, v] = convhulln(cube,'Pp'); +%! [h, v] = convhulln (cube); +%! assert (size (h), [6 4]); +%! h = sortrows (sort (h, 2), [1:4]); +%! assert (h, [1 2 3 4; 1 2 5 6; 1 4 5 8; 2 3 6 7; 3 4 7 8; 5 6 7 8]); +%! assert (v, 1, 10*eps); + +%!testif HAVE_QHULL +%! cube = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1]; +%! [h, v] = convhulln (cube, "QJ"); +%! assert (size (h), [12 3]); +%! assert (sortrows (sort (h, 2), [1:3]), [1 2 4; 1 2 5; 1 4 5; 2 3 4; 2 3 6; 2 5 6; 3 4 8; 3 6 7; 3 7 8; 4 5 8; 5 6 8; 6 7 8]); %! assert (v, 1.0, 1e6*eps); + %!testif HAVE_QHULL %! tetrahedron = [1 1 1;-1 -1 1;-1 1 -1;1 -1 -1]; -%! [h, v] = convhulln(tetrahedron,'Pp'); -%! assert (v, 8/3, 1e6*eps); +%! [h, v] = convhulln (tetrahedron); +%! h = sortrows (sort (h, 2), [1 2 3]); +%! assert (h, [1 2 3;1 2 4; 1 3 4; 2 3 4]); +%! assert (v, 8/3, 10*eps); */
--- a/src/DLD-FUNCTIONS/kron.cc Sat Oct 22 16:22:04 2011 -0700 +++ b/src/DLD-FUNCTIONS/kron.cc Tue Oct 25 11:28:41 2011 -0700 @@ -172,12 +172,77 @@ return octave_value (kron (am, bm)); } -#define ALL_TYPES(AMT, BMT) \ - } while (0) \ +octave_value +dispatch_kron (const octave_value& a, const octave_value& b) +{ + octave_value retval; + if (a.is_perm_matrix () && b.is_perm_matrix ()) + retval = do_kron<PermMatrix, PermMatrix> (a, b); + else if (a.is_diag_matrix ()) + { + if (b.is_diag_matrix () && a.rows () == a.columns () + && b.rows () == b.columns ()) + { + octave_value_list tmp; + tmp(0) = a.diag (); + tmp(1) = b.diag (); + tmp = dispatch_kron (tmp, 1); + if (tmp.length () == 1) + retval = tmp(0).diag (); + } + else if (a.is_single_type () || b.is_single_type ()) + { + if (a.is_complex_type ()) + retval = do_kron<FloatComplexDiagMatrix, FloatComplexMatrix> (a, b); + else if (b.is_complex_type ()) + retval = do_kron<FloatDiagMatrix, FloatComplexMatrix> (a, b); + else + retval = do_kron<FloatDiagMatrix, FloatMatrix> (a, b); + } + else + { + if (a.is_complex_type ()) + retval = do_kron<ComplexDiagMatrix, ComplexMatrix> (a, b); + else if (b.is_complex_type ()) + retval = do_kron<DiagMatrix, ComplexMatrix> (a, b); + else + retval = do_kron<DiagMatrix, Matrix> (a, b); + } + } + else if (a.is_sparse_type () || b.is_sparse_type ()) + { + if (a.is_complex_type () || b.is_complex_type ()) + retval = do_kron<SparseComplexMatrix, SparseComplexMatrix> (a, b); + else + retval = do_kron<SparseMatrix, SparseMatrix> (a, b); + } + else if (a.is_single_type () || b.is_single_type ()) + { + if (a.is_complex_type ()) + retval = do_kron<FloatComplexMatrix, FloatComplexMatrix> (a, b); + else if (b.is_complex_type ()) + retval = do_kron<FloatMatrix, FloatComplexMatrix> (a, b); + else + retval = do_kron<FloatMatrix, FloatMatrix> (a, b); + } + else + { + if (a.is_complex_type ()) + retval = do_kron<ComplexMatrix, ComplexMatrix> (a, b); + else if (b.is_complex_type ()) + retval = do_kron<Matrix, ComplexMatrix> (a, b); + else + retval = do_kron<Matrix, Matrix> (a, b); + } + return retval; +} + DEFUN_DLD (kron, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} kron (@var{A}, @var{B})\n\ -Form the Kronecker product of two matrices, defined block by block as\n\ +@deftypefnx {Loadable Function} {} kron (@var{A1}, @var{A2}, @dots{})\n\ +Form the Kronecker product of two or more matrices, defined block by \n\ +block as\n\ \n\ @example\n\ x = [a(i, j) b]\n\ @@ -193,86 +258,48 @@ 1 2 3 4\n\ @end group\n\ @end example\n\ +\n\ +If there are more than two input arguments @var{A1}, @var{A2}, @dots{}, \n\ +@var{An} the Kronecker product is computed as\n\ +\n\ +@example\n\ +kron (kron (@var{A1}, @var{A2}), @dots{}, @var{An})\n\ +@end example\n\ +\n\ +@noindent\n\ +Since the Kronecker product is associative, this is well-defined.\n\ @end deftypefn") { octave_value retval; int nargin = args.length (); - if (nargin == 2) + if (nargin >= 2) { octave_value a = args(0), b = args(1); - if (a.is_perm_matrix () && b.is_perm_matrix ()) - retval = do_kron<PermMatrix, PermMatrix> (a, b); - else if (a.is_diag_matrix ()) - { - if (b.is_diag_matrix () && a.rows () == a.columns () - && b.rows () == b.columns ()) - { - octave_value_list tmp; - tmp(0) = a.diag (); - tmp(1) = b.diag (); - tmp = Fkron (tmp, 1); - if (tmp.length () == 1) - retval = tmp(0).diag (); - } - else if (a.is_single_type () || b.is_single_type ()) - { - if (a.is_complex_type ()) - return do_kron<FloatComplexDiagMatrix, FloatComplexMatrix> (a, b); - else if (b.is_complex_type ()) - return do_kron<FloatDiagMatrix, FloatComplexMatrix> (a, b); - else - return do_kron<FloatDiagMatrix, FloatMatrix> (a, b); - } - else - { - if (a.is_complex_type ()) - return do_kron<ComplexDiagMatrix, ComplexMatrix> (a, b); - else if (b.is_complex_type ()) - return do_kron<DiagMatrix, ComplexMatrix> (a, b); - else - return do_kron<DiagMatrix, Matrix> (a, b); - } - } - else if (a.is_sparse_type () || b.is_sparse_type ()) - { - if (args(0).is_complex_type () || args(1).is_complex_type ()) - return do_kron<SparseComplexMatrix, SparseComplexMatrix> (a, b); - else - return do_kron<SparseMatrix, SparseMatrix> (a, b); - } - else if (a.is_single_type () || b.is_single_type ()) - { - if (a.is_complex_type ()) - return do_kron<FloatComplexMatrix, FloatComplexMatrix> (a, b); - else if (b.is_complex_type ()) - return do_kron<FloatMatrix, FloatComplexMatrix> (a, b); - else - return do_kron<FloatMatrix, FloatMatrix> (a, b); - } - else - { - if (a.is_complex_type ()) - return do_kron<ComplexMatrix, ComplexMatrix> (a, b); - else if (b.is_complex_type ()) - return do_kron<Matrix, ComplexMatrix> (a, b); - else - return do_kron<Matrix, Matrix> (a, b); - } + retval = dispatch_kron (a, b); + for (octave_idx_type i = 2; i < nargin; i++) + retval = dispatch_kron (retval, args(i)); } + else + print_usage (); return retval; } + /* %!test %! x = ones(2); %! assert( kron (x, x), ones (4)); -%!test +%!shared x, y, z +%! x = [1, 2]; +%! y = [-1, -2]; %! z = [1, 2, 3, 4; 1, 2, 3, 4; 1, 2, 3, 4]; -%! assert( kron (1:4, ones (3, 1)), z) +%!assert (kron (1:4, ones (3, 1)), z) +%!assert (kron (x, y, z), kron (kron (x, y), z)) +%!assert (kron (x, y, z), kron (x, kron (y, z))) */
--- a/src/help.cc Sat Oct 22 16:22:04 2011 -0700 +++ b/src/help.cc Tue Oct 25 11:28:41 2011 -0700 @@ -729,14 +729,14 @@ "-*- texinfo -*-\n\ @deftypefn {Keyword} {} varargin\n\ Pass an arbitrary number of arguments into a function.\n\ -@seealso{varargout, nargin, nargout}\n\ +@seealso{varargout, nargin, isargout, nargout, nthargout}\n\ @end deftypefn"), pair_type ("varargout", "-*- texinfo -*-\n\ @deftypefn {Keyword} {} varargout\n\ Pass an arbitrary number of arguments out of a function.\n\ -@seealso{varargin, nargin, nargout}\n\ +@seealso{varargin, nargin, isargout, nargout, nthargout}\n\ @end deftypefn"), pair_type ("while",
--- a/src/ov-usr-fcn.cc Sat Oct 22 16:22:04 2011 -0700 +++ b/src/ov-usr-fcn.cc Tue Oct 25 11:28:41 2011 -0700 @@ -635,7 +635,7 @@ Octave. If called with the optional argument @var{fcn_name}, return the\n\ maximum number of arguments the named function can accept, or -1 if the\n\ function accepts a variable number of arguments.\n\ -@seealso{nargout, varargin, varargout}\n\ +@seealso{nargout, varargin, isargout, varargout, nthargout}\n\ @end deftypefn") { octave_value retval; @@ -709,7 +709,7 @@ @code{f}.\n\ \n\ At the top level, @code{nargout} is undefined.\n\ -@seealso{nargin, isargout, varargin, varargout}\n\ +@seealso{nargin, varargin, isargout, varargout, nthargout}\n\ @end deftypefn") { octave_value retval; @@ -806,7 +806,7 @@ false. @var{k} can also be an array, in which case the function works\n\ element-by-element and a logical array is returned. At the top level,\n\ @code{isargout} returns an error.\n\ -@seealso{nargout, nargin, varargin, varargout}\n\ +@seealso{nargout, nargin, varargin, varargout, nthargout}\n\ @end deftypefn") { octave_value retval;