Mercurial > octave-nkf
changeset 9899:9f25290a35e8
more private function and subfunction changes
line wrap: on
line diff
--- a/scripts/ChangeLog Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/ChangeLog Tue Dec 01 22:40:37 2009 -0500 @@ -1,3 +1,59 @@ +2009-12-01 John W. Eaton <jwe@octave.org> + + * help/module.mk (help_PRIVATE_FCN_FILES): New list. + (help_FCN_FILES): Remove new private functions from the list. + Include $(help_PRIVATE_FCN_FILES) in the list. + * help/private/__additional_help_message__.m: Rename from + help/__additional_help_message__.m. + + * statistics/base/module.mk (statistics_base_FCN_FILES): + Remove statistics/base/__quantile__.m from the list. + * statistics/base/__quantile__.m: Now a subfunction of + statistics/base/quantile.m. + * statistics/base/quantile.m: Remove redundant tests. + + * miscellaneous/__xzip__.m: Comment out tests until we have a way + to test private functions directly. + + * general/isequal.m, general/isequalwithequalnans.m: + Convert tests from __isequal__. + + * optimization/module.mk (optimization_PRIVATE_FCN_FILES): New list. + (optimization_FCN_FILES): Remove new private functions and new + subfunctions from the list. Include + $(optimization_PRIVATE_FCN_FILES) in the list. + + * optimization/private/__fdjac__.m: Rename from + optimization/__fdjac__.m. + + * optimization/__dogleg__.m: Now a subfunction of path/fsolve.m. + * optimization/__doglegm__.m: Now a subfunction of path/fminunc.m. + + * general/module.mk (general_PRIVATE_FCN_FILES): New list. + (general_FCN_FILES): Remove new private functions from the list. + Include $(general_PRIVATE_FCN_FILES) in the list. + + * general/private/__isequal__.m: Rename from general/__isequal__.m. + * general/private/__splinen__.m: Rename from general/__splinen__.m. + + * image/module.mk (image_FCN_FILES): Remove image/__img__.m and + image/__img_via_file__.m from the list. + + * image/__img__.m: Now a subfunction of image/image.m. + * image/__img_via_file__.m: Now a subfunction of image_viewer.m. + + * path/module.mk (path_FCN_FILES): Remove path/__extractpath__.m + from the list. + + * path/__extractpath__.m: Now a subfunction of path/pathdef.m. + + * miscellaneous/module.mk (miscellaneous_PRIVATE_FCN_FILES): New list. + (miscellaneous_FCN_FILES): Remove __xzip__.m from the list. + Include $(miscellaneous_PRIVATE_FCN_FILES) in the list. + + * miscellaneous/private/__xzip__.m: Rename from + miscellaneous/__xzip__.m. + 2009-12-01 David Bateman <dbateman@free.fr> * @ftp/ftp.m: Treat empty constructor and construction from @@ -8,6 +64,10 @@ 2009-12-01 John W. Eaton <jwe@octave.org> + * plot/module.mk (plot_PRIVATE_FCN_FILES): New list. + (plot_FCN_FILES): Include $(plot_PRIVATE_FCN_FILES) in the list. + Remove new private functions and new subfunctions from the list. + * plot/private/__actual_axis_position__.m: Rename from plot/__actual_axis_position__.m. * plot/private/__add_datasource__.m: Rename from
--- a/scripts/general/__isequal__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,230 +0,0 @@ -## Copyright (C) 2000, 2005, 2006, 2007, 2009 Paul Kienzle -## -## 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/>. - -## Undocumented internal function. - -## -*- texinfo -*- -## @deftypefn {Function File} {} __isequal__ (@var{nans_compare_equal}, @var{x1}, @var{x2}, @dots{}) -## Undocumented internal function. -## @end deftypefn - -## Return true if @var{x1}, @var{x2}, @dots{} are all equal and -## @var{nans_compare_equal} evaluates to false. -## -## If @var{nans_compare_equal} evaluates to true, then assume NaN == NaN. - -## Modified by: William Poetra Yoga Hadisoeseno - -## Algorithm: -## -## 1. Determine the class of x -## 2. If x is of the struct, cell, list or char class, for each -## argument after x, determine whether it has the same class -## and size as x. -## Otherwise, for each argument after x, verify that it is not -## of the struct, cell, list or char class, and that it has -## the same size as x. -## 3. For each argument after x, compare it for equality with x: -## a. struct compare each member by name, not by order (recursive) -## b. cell/list compare each member by order (recursive) -## c. char compare each member with strcmp -## d. <other> compare each nonzero member, and assume NaN == NaN -## if nans_compare_equal is nonzero. - -function t = __isequal__ (nans_compare_equal, x, varargin) - - if (nargin < 3) - print_usage (); - endif - - l_v = nargin - 2; - - ## Generic tests. - - ## Give an error for a list (that will make the code simpler and lists - ## are deprecated anyway. - if (islist (x)) - error ("__isequal__: list objects are deprecated and cannot be tested for equality here; use cell arrays instead"); - endif - - ## All arguments must either be of the same class or they must be - ## numeric values. - t = (all (strcmp (class(x), - cellfun (@class, varargin, "UniformOutput", false))) - || ((isnumeric (x) || islogical (x)) - && all ((cellfun (@isnumeric, varargin) | cellfun (@islogical, varargin))))); - - if (t) - ## Test that everything has the same number of dimensions. - s_x = size (x); - s_v = cellfun (@size, varargin, "UniformOutput", false); - t = all (length (s_x) == cellfun (@length, s_v)); - endif - - if (t) - ## Test that everything is the same size since it has the same - ## dimensionality. - l_x = length (s_x); - s_v = reshape ([s_v{:}], length (s_x), []); - idx = 0; - while (t && idx < l_x) - idx++; - t = all (s_x(idx) == s_v(idx, :)); - endwhile - endif - - if (t) - ## Check individual classes. - if (isstruct (x)) - ## Test the number of fields. - fn_x = fieldnames (x); - l_fn_x = length (fn_x); - fn_v = cellfun (@fieldnames, varargin, "UniformOutput", false); - t = all (l_fn_x == cellfun (@length, fn_v)); - - ## Test that all the names are equal. - idx = 0; - s_fn_x = sort (fn_x); - while (t && idx < l_v) - idx++; - ## We'll allow the fieldnames to be in a different order. - t = all (strcmp (s_fn_x, sort (fn_v{idx}))); - endwhile - - idx = 0; - while (t && idx < l_fn_x) - ## Test that all field values are equal. - idx++; - args = {nans_compare_equal, {x.(fn_x{idx})}}; - for argn = 1:l_v - args{argn+2} = {varargin{argn}.(fn_x{idx})}; - endfor - ## Minimize function calls by calling for all the arguments at - ## once. - t = __isequal__ (args{:}); - endwhile - - elseif (iscell (x)) - ## Check that each element of a cell is equal. - l_x = numel (x); - idx = 0; - while (t && idx < l_x) - idx++; - args = {nans_compare_equal, x{idx}}; - for p = 1:l_v - args{p+2} = varargin{p}{idx}; - endfor - t = __isequal__ (args{:}); - endwhile - - elseif (ischar (x)) - - ## Sizes are equal already, so we can just make everything into a - ## row and test the rows. - for i = 1:l_v - strings{i} = reshape (varargin{i}, 1, []); - endfor - t = all (strcmp (reshape (x, 1, []), strings)); - - else - ## Check the numeric types. - - if (issparse (x)) - f_x = spfind (x); - else - f_x = find (x); - endif - l_f_x = length (f_x); - x = x(f_x); - for argn = 1:l_v - y = varargin{argn}; - if (issparse (y)) - f_y = spfind (y); - else - f_y = find (y); - endif - - t = (l_f_x == length (f_y)) && all (f_x == f_y); - if (!t) - return; - endif - - y = y(f_y); - m = (x == y); - t = all (m); - - if (!t) - if (nans_compare_equal) - t = isnan (x(!m)) && isnan (y(!m)); - else - return; - endif - endif - endfor - - endif - endif - -endfunction - -## test size and shape -%!assert(__isequal__(0,[1,2,3,4],[1,2,3,4]), true) -%!assert(__isequal__(0,[1;2;3;4],[1;2;3;4]), true) -%!assert(__isequal__(0,[1,2,3,4],[1;2;3;4]), false) -%!assert(__isequal__(0,[1,2,3,4],[1,2;3,4]), false) -%!assert(__isequal__(0,[1,2,3,4],[1,3;2,4]), false) - -%!test -%! A = 1:8; -%! B = reshape (A, 2, 2, 2); -%! assert (__isequal__ (0, A, B), false); - -%!test -%! A = reshape (1:8, 2, 2, 2); -%! B = A; -%! assert (__isequal__ (0, A, B), true); - -%!test -%! A = reshape (1:8, 2, 4); -%! B = reshape (A, 2, 2, 2); -%! assert (__isequal__ (0, A, B), false); - -## test for equality -%!assert(__isequal__(0,[1,2,3,4],[1,2,3,4]), true) -%!assert(__isequal__(1,{1,2,NaN,4},{1,2,NaN,4}), true) -%!assert(__isequal__(1,[1,2,NaN,4],[1,2,NaN,4]), true) -%!assert(__isequal__(0,['a','b','c','d'],['a','b','c','d']), true) -## Test multi-line strings -%!assert(__isequal__(0,["test";"strings"],["test";"strings"],["test";"strings"]), true) -## test for inequality -%!assert(__isequal__(0,[1,2,3,4],[1;2;3;4]),false) -%!assert(__isequal__(0,{1,2,3,4},[1,2,3,4]),false) -%!assert(__isequal__(0,[1,2,3,4],{1,2,3,4}),false) -%!assert(__isequal__(0,[1,2,NaN,4],[1,2,NaN,4]),false) -%!assert(__isequal__(1,[1,2,NaN,4],[1,NaN,3,4]),false) -%!assert(__isequal__(1,[1,2,NaN,4],[1,2,3,4]),false) -%!assert(__isequal__(0,['a','b','c','d'],['a';'b';'c';'d']),false) -%!assert(__isequal__(0,{'a','b','c','d'},{'a';'b';'c';'d'}),false) -## test for equality (struct) -%!assert(__isequal__(0,struct('a',1,'b',2),struct('a',1,'b',2)),true) -%!assert(__isequal__(0,struct('a',1,'b',2),struct('a',1,'b',2),struct('a',1,'b',2)),true) -%!assert(__isequal__(0,struct('a','abc','b',2),struct('a','abc','b',2)),true) -%!assert(__isequal__(1,struct('a',NaN,'b',2),struct('a',NaN,'b',2),struct('a',NaN,'b',2)),true) -## test for inequality (struct) -%!assert(__isequal__(0,struct('a',NaN,'b',2),struct('a',NaN,'b',2),struct('a',NaN,'b',2)),false) -
--- a/scripts/general/__splinen__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,49 +0,0 @@ -## Copyright (C) 2007, 2008, 2009 David Bateman -## -## 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/>. - -## Undocumented internal function. - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{yi} =} __splinen__ (@var{x}, @var{y}, @var{xi}) -## Undocumented internal function. -## @end deftypefn - -## FIXME: Allow arbitrary grids.. - -function yi = __splinen__ (x, y, xi, extrapval, f) - if (nargin != 5) - error ("Incorrect number of arguments"); - endif - ## ND isvector function. - isvec = @(x) numel (x) == length (x); - if (!iscell (x) || length(x) < ndims(y) || any (! cellfun (isvec, x)) || - !iscell (xi) || length(xi) < ndims(y) || any (! cellfun (isvec, xi))) - error ("%s: non gridded data or dimensions inconsistent", f); - endif - yi = y; - for i = length(x):-1:1 - yi = permute (spline (x{i}, yi, xi{i}(:)), [length(x),1:length(x)-1]); - endfor - - [xi{:}] = ndgrid (cellfun (@(x) x(:), xi, "UniformOutput", false){:}); - idx = zeros (size(xi{1})); - for i = 1 : length(x) - idx |= xi{i} < min (x{i}(:)) | xi{i} > max (x{i}(:)); - endfor - yi(idx) = extrapval; -endfunction
--- a/scripts/general/isequal.m Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/general/isequal.m Tue Dec 01 22:40:37 2009 -0500 @@ -25,10 +25,50 @@ function retval = isequal (x, varargin) if (nargin > 1) - retval = __isequal__ (0, x, varargin{:}); + retval = __isequal__ (false, x, varargin{:}); else print_usage (); endif endfunction +## test size and shape +%!assert(isequal([1,2,3,4],[1,2,3,4]), true) +%!assert(isequal([1;2;3;4],[1;2;3;4]), true) +%!assert(isequal([1,2,3,4],[1;2;3;4]), false) +%!assert(isequal([1,2,3,4],[1,2;3,4]), false) +%!assert(isequal([1,2,3,4],[1,3;2,4]), false) + +%!test +%! A = 1:8; +%! B = reshape (A, 2, 2, 2); +%! assert (isequal (A, B), false); + +%!test +%! A = reshape (1:8, 2, 2, 2); +%! B = A; +%! assert (isequal (A, B), true); + +%!test +%! A = reshape (1:8, 2, 4); +%! B = reshape (A, 2, 2, 2); +%! assert (isequal (A, B), false); + +## test for equality +%!assert(isequal([1,2,3,4],[1,2,3,4]), true) +%!assert(isequal(['a','b','c','d'],['a','b','c','d']), true) +## Test multi-line strings +%!assert(isequal(["test";"strings"],["test";"strings"],["test";"strings"]), true) +## test for inequality +%!assert(isequal([1,2,3,4],[1;2;3;4]),false) +%!assert(isequal({1,2,3,4},[1,2,3,4]),false) +%!assert(isequal([1,2,3,4],{1,2,3,4}),false) +%!assert(isequal([1,2,NaN,4],[1,2,NaN,4]),false) +%!assert(isequal(['a','b','c','d'],['a';'b';'c';'d']),false) +%!assert(isequal({'a','b','c','d'},{'a';'b';'c';'d'}),false) +## test for equality (struct) +%!assert(isequal(struct('a',1,'b',2),struct('a',1,'b',2)),true) +%!assert(isequal(struct('a',1,'b',2),struct('a',1,'b',2),struct('a',1,'b',2)),true) +%!assert(isequal(struct('a','abc','b',2),struct('a','abc','b',2)),true) +## test for inequality (struct) +%!assert(isequal(struct('a',NaN,'b',2),struct('a',NaN,'b',2),struct('a',NaN,'b',2)),false)
--- a/scripts/general/isequalwithequalnans.m Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/general/isequalwithequalnans.m Tue Dec 01 22:40:37 2009 -0500 @@ -26,10 +26,18 @@ function retval = isequalwithequalnans (x, varargin) if (nargin > 1) - retval = __isequal__ (1, x, varargin{:}); + retval = __isequal__ (true, x, varargin{:}); else print_usage (); endif endfunction +## test for equality +%!assert(isequalwithequalnans({1,2,NaN,4},{1,2,NaN,4}), true) +%!assert(isequalwithequalnans([1,2,NaN,4],[1,2,NaN,4]), true) +## test for inequality +%!assert(isequalwithequalnans([1,2,NaN,4],[1,NaN,3,4]),false) +%!assert(isequalwithequalnans([1,2,NaN,4],[1,2,3,4]),false) +## test for equality (struct) +%!assert(isequalwithequalnans(struct('a',NaN,'b',2),struct('a',NaN,'b',2),struct('a',NaN,'b',2)),true)
--- a/scripts/general/module.mk Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/general/module.mk Tue Dec 01 22:40:37 2009 -0500 @@ -1,8 +1,10 @@ FCN_FILE_DIRS += general +general_PRIVATE_FCN_FILES = \ + general/private/__isequal__.m \ + general/private/__splinen__.m + general_FCN_FILES = \ - general/__isequal__.m \ - general/__splinen__.m \ general/accumarray.m \ general/arrayfun.m \ general/bicubic.m \ @@ -75,7 +77,8 @@ general/structfun.m \ general/subsindex.m \ general/triplequad.m \ - general/trapz.m + general/trapz.m \ + $(general_PRIVATE_FCN_FILES) FCN_FILES += $(general_FCN_FILES)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/private/__isequal__.m Tue Dec 01 22:40:37 2009 -0500 @@ -0,0 +1,183 @@ +## Copyright (C) 2000, 2005, 2006, 2007, 2009 Paul Kienzle +## +## 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/>. + +## Undocumented internal function. + +## -*- texinfo -*- +## @deftypefn {Function File} {} __isequal__ (@var{nans_compare_equal}, @var{x1}, @var{x2}, @dots{}) +## Undocumented internal function. +## @end deftypefn + +## Return true if @var{x1}, @var{x2}, @dots{} are all equal and +## @var{nans_compare_equal} evaluates to false. +## +## If @var{nans_compare_equal} evaluates to true, then assume NaN == NaN. + +## Modified by: William Poetra Yoga Hadisoeseno + +## Algorithm: +## +## 1. Determine the class of x +## 2. If x is of the struct, cell, list or char class, for each +## argument after x, determine whether it has the same class +## and size as x. +## Otherwise, for each argument after x, verify that it is not +## of the struct, cell, list or char class, and that it has +## the same size as x. +## 3. For each argument after x, compare it for equality with x: +## a. struct compare each member by name, not by order (recursive) +## b. cell/list compare each member by order (recursive) +## c. char compare each member with strcmp +## d. <other> compare each nonzero member, and assume NaN == NaN +## if nans_compare_equal is nonzero. + +function t = __isequal__ (nans_compare_equal, x, varargin) + + if (nargin < 3) + print_usage (); + endif + + l_v = nargin - 2; + + ## Generic tests. + + ## Give an error for a list (that will make the code simpler and lists + ## are deprecated anyway. + if (islist (x)) + error ("__isequal__: list objects are deprecated and cannot be tested for equality here; use cell arrays instead"); + endif + + ## All arguments must either be of the same class or they must be + ## numeric values. + t = (all (strcmp (class(x), + cellfun (@class, varargin, "UniformOutput", false))) + || ((isnumeric (x) || islogical (x)) + && all ((cellfun (@isnumeric, varargin) | cellfun (@islogical, varargin))))); + + if (t) + ## Test that everything has the same number of dimensions. + s_x = size (x); + s_v = cellfun (@size, varargin, "UniformOutput", false); + t = all (length (s_x) == cellfun (@length, s_v)); + endif + + if (t) + ## Test that everything is the same size since it has the same + ## dimensionality. + l_x = length (s_x); + s_v = reshape ([s_v{:}], length (s_x), []); + idx = 0; + while (t && idx < l_x) + idx++; + t = all (s_x(idx) == s_v(idx, :)); + endwhile + endif + + if (t) + ## Check individual classes. + if (isstruct (x)) + ## Test the number of fields. + fn_x = fieldnames (x); + l_fn_x = length (fn_x); + fn_v = cellfun (@fieldnames, varargin, "UniformOutput", false); + t = all (l_fn_x == cellfun (@length, fn_v)); + + ## Test that all the names are equal. + idx = 0; + s_fn_x = sort (fn_x); + while (t && idx < l_v) + idx++; + ## We'll allow the fieldnames to be in a different order. + t = all (strcmp (s_fn_x, sort (fn_v{idx}))); + endwhile + + idx = 0; + while (t && idx < l_fn_x) + ## Test that all field values are equal. + idx++; + args = {nans_compare_equal, {x.(fn_x{idx})}}; + for argn = 1:l_v + args{argn+2} = {varargin{argn}.(fn_x{idx})}; + endfor + ## Minimize function calls by calling for all the arguments at + ## once. + t = __isequal__ (args{:}); + endwhile + + elseif (iscell (x)) + ## Check that each element of a cell is equal. + l_x = numel (x); + idx = 0; + while (t && idx < l_x) + idx++; + args = {nans_compare_equal, x{idx}}; + for p = 1:l_v + args{p+2} = varargin{p}{idx}; + endfor + t = __isequal__ (args{:}); + endwhile + + elseif (ischar (x)) + + ## Sizes are equal already, so we can just make everything into a + ## row and test the rows. + for i = 1:l_v + strings{i} = reshape (varargin{i}, 1, []); + endfor + t = all (strcmp (reshape (x, 1, []), strings)); + + else + ## Check the numeric types. + + if (issparse (x)) + f_x = spfind (x); + else + f_x = find (x); + endif + l_f_x = length (f_x); + x = x(f_x); + for argn = 1:l_v + y = varargin{argn}; + if (issparse (y)) + f_y = spfind (y); + else + f_y = find (y); + endif + + t = (l_f_x == length (f_y)) && all (f_x == f_y); + if (!t) + return; + endif + + y = y(f_y); + m = (x == y); + t = all (m); + + if (!t) + if (nans_compare_equal) + t = isnan (x(!m)) && isnan (y(!m)); + else + return; + endif + endif + endfor + + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/private/__splinen__.m Tue Dec 01 22:40:37 2009 -0500 @@ -0,0 +1,49 @@ +## Copyright (C) 2007, 2008, 2009 David Bateman +## +## 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/>. + +## Undocumented internal function. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{yi} =} __splinen__ (@var{x}, @var{y}, @var{xi}) +## Undocumented internal function. +## @end deftypefn + +## FIXME: Allow arbitrary grids.. + +function yi = __splinen__ (x, y, xi, extrapval, f) + if (nargin != 5) + error ("Incorrect number of arguments"); + endif + ## ND isvector function. + isvec = @(x) numel (x) == length (x); + if (!iscell (x) || length(x) < ndims(y) || any (! cellfun (isvec, x)) || + !iscell (xi) || length(xi) < ndims(y) || any (! cellfun (isvec, xi))) + error ("%s: non gridded data or dimensions inconsistent", f); + endif + yi = y; + for i = length(x):-1:1 + yi = permute (spline (x{i}, yi, xi{i}(:)), [length(x),1:length(x)-1]); + endfor + + [xi{:}] = ndgrid (cellfun (@(x) x(:), xi, "UniformOutput", false){:}); + idx = zeros (size(xi{1})); + for i = 1 : length(x) + idx |= xi{i} < min (x{i}(:)) | xi{i} > max (x{i}(:)); + endfor + yi(idx) = extrapval; +endfunction
--- a/scripts/help/__additional_help_message__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,37 +0,0 @@ -## Copyright (C) 2009 Søren Hauberg -## -## This program 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. -## -## This program 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 this program; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn {Function File} {} __additional_help_message__ () -## Undocumented internal function. -## @end deftypefn - -function msg = __additional_help_message__ () - - if (suppress_verbose_help_message ()) - msg = ""; - else - msg = "\ -Additional help for built-in functions and operators is\n\ -available in the on-line version of the manual. Use the command\n\ -`doc <topic>' to search the manual index.\n\ -\n\ -Help and information about Octave is also available on the WWW\n\ -at http://www.octave.org and via the help@octave.org\n\ -mailing list.\n"; - endif - -endfunction
--- a/scripts/help/module.mk Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/help/module.mk Tue Dec 01 22:40:37 2009 -0500 @@ -1,7 +1,9 @@ FCN_FILE_DIRS += help +help_PRIVATE_FCN_FILES = \ + help/__additional_help_message__.m \ + help_FCN_FILES = \ - help/__additional_help_message__.m \ help/__makeinfo__.m \ help/__strip_html_tags__.m \ help/doc.m \ @@ -11,7 +13,8 @@ help/lookfor.m \ help/print_usage.m \ help/type.m \ - help/which.m + help/which.m \ + $(help_PRIVATE_FCN_FILES) FCN_FILES += $(help_FCN_FILES)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/help/private/__additional_help_message__.m Tue Dec 01 22:40:37 2009 -0500 @@ -0,0 +1,37 @@ +## Copyright (C) 2009 Søren Hauberg +## +## This program 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. +## +## This program 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 this program; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} __additional_help_message__ () +## Undocumented internal function. +## @end deftypefn + +function msg = __additional_help_message__ () + + if (suppress_verbose_help_message ()) + msg = ""; + else + msg = "\ +Additional help for built-in functions and operators is\n\ +available in the on-line version of the manual. Use the command\n\ +`doc <topic>' to search the manual index.\n\ +\n\ +Help and information about Octave is also available on the WWW\n\ +at http://www.octave.org and via the help@octave.org\n\ +mailing list.\n"; + endif + +endfunction
--- a/scripts/image/__img__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,90 +0,0 @@ -## Copyright (C) 1996, 1997, 2006, 2007, 2008, 2009 John W. Eaton -## -## 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/>. - -## Undocumented internal function. - -## -*- texinfo -*- -## @deftypefnx {Function File} {} __img__ (@var{x}, @var{y}, @var{img}, @dots{}) -## Undocumented internal function. -## @end deftypefn - -## Generic image creation. -## -## The axis values corresponding to the matrix elements are specified in -## @var{x} and @var{y}. If you're not using gnuplot 4.2 or later, these -## variables are ignored. - -## Author: Tony Richardson <arichard@stark.cc.oh.us> -## Created: July 1994 -## Adapted-By: jwe - -function h = __img__ (x, y, img, varargin) - - newplot (); - - if (isempty (img)) - error ("__img__: matrix is empty"); - endif - - if (isempty (x)) - x = [1, columns(img)]; - endif - - if (isempty (y)) - y = [1, rows(img)]; - endif - - xdata = [x(1), x(end)]; - ydata = [y(1), y(end)]; - - xlim = [x(1)-0.5, x(end)+0.5]; - ylim = [y(1)-0.5, y(end)+0.5]; - - ca = gca (); - - tmp = __go_image__ (ca, "cdata", img, "xdata", xdata, "ydata", ydata, - "cdatamapping", "direct", varargin {:}); - - ## FIXME -- how can we do this and also get the {x,y}limmode - ## properties to remain "auto"? I suppose this adjustment should - ## happen automatically in axes::update_axis_limits instead of - ## explicitly setting the values here. But then what information is - ## available to axes::update_axis_limits to determine that the - ## adjustment is necessary? - set (ca, "xlim", xlim, "ylim", ylim); - - if (ndims (img) == 3) - if (isinteger (img)) - c = class (img); - mn = intmin (c); - mx = intmax (c); - set (ca, "clim", double ([mn, mx])); - endif - endif - - set (ca, "view", [0, 90]); - - if (strcmp (get (ca, "nextplot"), "replace")) - set (ca, "ydir", "reverse"); - endif - - if (nargout > 0) - h = tmp; - endif - -endfunction
--- a/scripts/image/__img_via_file__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,66 +0,0 @@ -## Copyright (C) 2006, 2007, 2009 Søren Hauberg -## -## 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/>. - -## Undocumented internal function. - -## -*- texinfo -*- -## @deftypefn {Function File} {} __img_via_file__ (@var{x}, @var{y}, @var{im}, @var{zoom}, @var{command}) -## Undocumented internal function. -## @end deftypefn - -## Display an image by saving it to a file in PPM format and launching -## @var{command}. -## -## The @var{command} must be a format string containing @code{%s} and -## possibly @code{%f}. The @code{%s} will be replaced by the filename -## of the image, and the @code{%f} will be replaced by @var{zoom}. The -## @var{x} and @var{y} arguments are ignored. - -function __img_via_file__ (x, y, im, zoom, command) - - ppm_name = tmpnam (); - saveimage (ppm_name, im, "ppm"); - - rm = sprintf ("rm -f \"%s\"", ppm_name); - - if (isempty (command)) - ## Different image viewer commands. - xv = sprintf ("xv -raw -expand %f \"%s\"", zoom, ppm_name); - xloadimage = sprintf ("xloadimage -zoom %f \"%s\"", zoom*100, ppm_name); - im_display = sprintf ("display -resize %f%% \"%s\"", zoom*100, ppm_name); - - ## Need to let the shell clean up the tmp file because we are putting - ## the viewer in the background. - status = system (sprintf ("( %s || %s || %s && %s ) > /dev/null 2>&1 &", - im_display, xv, xloadimage, rm)); - else - ## Does the command support zooming? - if (findstr (command, "%f")) - command = sprintf (command, zoom, ppm_name); - else - command = sprintf (command, ppm_name); - endif - status = system (sprintf ("( %s && %s) > /dev/null 2>&1 &", command, rm)); - endif - - ## Did the system call fail? - if (status != 0) - error ("the image viewing command failed"); - endif - -endfunction
--- a/scripts/image/image.m Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/image/image.m Tue Dec 01 22:40:37 2009 -0500 @@ -80,3 +80,69 @@ endif endfunction + +## Generic image creation. +## +## The axis values corresponding to the matrix elements are specified in +## @var{x} and @var{y}. If you're not using gnuplot 4.2 or later, these +## variables are ignored. + +## Author: Tony Richardson <arichard@stark.cc.oh.us> +## Created: July 1994 +## Adapted-By: jwe + +function h = __img__ (x, y, img, varargin) + + newplot (); + + if (isempty (img)) + error ("__img__: matrix is empty"); + endif + + if (isempty (x)) + x = [1, columns(img)]; + endif + + if (isempty (y)) + y = [1, rows(img)]; + endif + + xdata = [x(1), x(end)]; + ydata = [y(1), y(end)]; + + xlim = [x(1)-0.5, x(end)+0.5]; + ylim = [y(1)-0.5, y(end)+0.5]; + + ca = gca (); + + tmp = __go_image__ (ca, "cdata", img, "xdata", xdata, "ydata", ydata, + "cdatamapping", "direct", varargin {:}); + + ## FIXME -- how can we do this and also get the {x,y}limmode + ## properties to remain "auto"? I suppose this adjustment should + ## happen automatically in axes::update_axis_limits instead of + ## explicitly setting the values here. But then what information is + ## available to axes::update_axis_limits to determine that the + ## adjustment is necessary? + set (ca, "xlim", xlim, "ylim", ylim); + + if (ndims (img) == 3) + if (isinteger (img)) + c = class (img); + mn = intmin (c); + mx = intmax (c); + set (ca, "clim", double ([mn, mx])); + endif + endif + + set (ca, "view", [0, 90]); + + if (strcmp (get (ca, "nextplot"), "replace")) + set (ca, "ydir", "reverse"); + endif + + if (nargout > 0) + h = tmp; + endif + +endfunction
--- a/scripts/image/image_viewer.m Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/image/image_viewer.m Tue Dec 01 22:40:37 2009 -0500 @@ -119,3 +119,45 @@ endif endfunction + +## Display an image by saving it to a file in PPM format and launching +## @var{command}. +## +## The @var{command} must be a format string containing @code{%s} and +## possibly @code{%f}. The @code{%s} will be replaced by the filename +## of the image, and the @code{%f} will be replaced by @var{zoom}. The +## @var{x} and @var{y} arguments are ignored. + +function __img_via_file__ (x, y, im, zoom, command) + + ppm_name = tmpnam (); + saveimage (ppm_name, im, "ppm"); + + rm = sprintf ("rm -f \"%s\"", ppm_name); + + if (isempty (command)) + ## Different image viewer commands. + xv = sprintf ("xv -raw -expand %f \"%s\"", zoom, ppm_name); + xloadimage = sprintf ("xloadimage -zoom %f \"%s\"", zoom*100, ppm_name); + im_display = sprintf ("display -resize %f%% \"%s\"", zoom*100, ppm_name); + + ## Need to let the shell clean up the tmp file because we are putting + ## the viewer in the background. + status = system (sprintf ("( %s || %s || %s && %s ) > /dev/null 2>&1 &", + im_display, xv, xloadimage, rm)); + else + ## Does the command support zooming? + if (findstr (command, "%f")) + command = sprintf (command, zoom, ppm_name); + else + command = sprintf (command, ppm_name); + endif + status = system (sprintf ("( %s && %s) > /dev/null 2>&1 &", command, rm)); + endif + + ## Did the system call fail? + if (status != 0) + error ("the image viewing command failed"); + endif + +endfunction
--- a/scripts/image/module.mk Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/image/module.mk Tue Dec 01 22:40:37 2009 -0500 @@ -1,8 +1,6 @@ FCN_FILE_DIRS += image image_FCN_FILES = \ - image/__img__.m \ - image/__img_via_file__.m \ image/autumn.m \ image/bone.m \ image/brighten.m \
--- a/scripts/miscellaneous/__xzip__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,140 +0,0 @@ -## Copyright (C) 2008, 2009 Thorsten Meyer -## based on the original gzip function by David Bateman -## -## 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} {@var{entries} =} __xzip__ (@var{commandname}, @var{extension}, @var{commandtemplate}, @var{files}, @var{outdir}) -## Undocumented internal function. -## @end deftypefn - -## Compress the list of files and/or directories specified in @var{files} -## with the external compression command @var{commandname}. The template -## @var{commandtemplate} is used to actually start the command. Each file -## is compressed separately and a new file with the extension @var{extension} -## is created and placed into the directory @var{outdir}. The original files -## are not touched. Existing compressed files are silently overwritten. -## This is an internal function. Do not use directly. - -function entries = __xzip__ (commandname, extension, - commandtemplate, files, outdir) - - if (nargin == 4 || nargin == 5) - if (! ischar (extension) || length (extension) == 0) - error ("__xzip__: extension has to be a string with finite length"); - endif - - if (nargin == 5 && ! exist (outdir, "dir")) - error ("__xzip__: output directory does not exist"); - endif - - if (ischar (files)) - files = cellstr (files); - else - error ("__xzip__: expecting FILES to be a character array"); - endif - - if (nargin == 4) - outdir = tmpnam (); - mkdir (outdir); - endif - - cwd = pwd(); - unwind_protect - files = glob (files); - - ## Ignore any file with the compress extension - files (cellfun (@(x) length(x) > length(extension) - && strcmp (x((end - length(extension) + 1):end), extension), - files)) = []; - - copyfile (files, outdir); - - [d, f] = myfileparts(files); - - cd (outdir); - - cmd = sprintf (commandtemplate, sprintf (" %s", f{:})); - - [status, output] = system (cmd); - if (status == 0) - - if (nargin == 5) - compressed_files = cellfun( - @(x) fullfile (outdir, sprintf ("%s.%s", x, extension)), - f, "UniformOutput", false); - else - movefile (cellfun(@(x) sprintf ("%s.%s", x, extension), f, - "UniformOutput", false), cwd); - ## FIXME this does not work when you try to compress directories - - compressed_files = cellfun(@(x) sprintf ("%s.%s", x, extension), - files, "UniformOutput", false); - endif - - if (nargout > 0) - entries = compressed_files; - endif - else - error ("%s command failed with exit status = %d", - commandname, status); - endif - unwind_protect_cleanup - cd(cwd); - if (nargin == 4) - crr = confirm_recursive_rmdir (); - unwind_protect - confirm_recursive_rmdir (false); - rmdir (outdir, "s"); - unwind_protect_cleanup - confirm_recursive_rmdir (crr); - end_unwind_protect - endif - end_unwind_protect - else - print_usage (); - endif - -endfunction - -function [d, f] = myfileparts (files) - [d, f, ext] = cellfun (@(x) fileparts (x), files, "UniformOutput", false); - f = cellfun (@(x, y) sprintf ("%s%s", x, y), f, ext, - "UniformOutput", false); - idx = cellfun (@(x) isdir (x), files); - d(idx) = ""; - f(idx) = files(idx); -endfunction - -%!error <extension has to be a string with finite length> -%! __xzip__("gzip", "", "gzip -r %s", "bla"); -%!error <no files to move> -%! __xzip__("gzip", ".gz", "gzip -r %s", tmpnam); -%!error <command failed with exit status> -%! # test __xzip__ with invalid compression command -%! unwind_protect -%! filename = tmpnam; -%! dummy = 1; -%! save(filename, "dummy"); -%! dirname = tmpnam; -%! mkdir(dirname); -%! entry = __xzip__("gzip", ".gz", "xxxzipxxx -r %s 2>/dev/null", -%! filename, dirname); -%! unwind_protect_cleanup -%! delete(filename); -%! rmdir(dirname); -%! end_unwind_protect
--- a/scripts/miscellaneous/module.mk Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/miscellaneous/module.mk Tue Dec 01 22:40:37 2009 -0500 @@ -1,7 +1,9 @@ FCN_FILE_DIRS += miscellaneous +miscellaneous_PRIVATE_FCN_FILES = \ + miscellaneous/private/__xzip__.m + miscellaneous_FCN_FILES = \ - miscellaneous/__xzip__.m \ miscellaneous/ans.m \ miscellaneous/bincoeff.m \ miscellaneous/bug_report.m \ @@ -66,7 +68,8 @@ miscellaneous/warning_ids.m \ miscellaneous/what.m \ miscellaneous/xor.m \ - miscellaneous/zip.m + miscellaneous/zip.m \ + $(miscellaneous_PRIVATE_FCN_FILES) FCN_FILES += $(miscellaneous_FCN_FILES)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/miscellaneous/private/__xzip__.m Tue Dec 01 22:40:37 2009 -0500 @@ -0,0 +1,143 @@ +## Copyright (C) 2008, 2009 Thorsten Meyer +## based on the original gzip function by David Bateman +## +## 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} {@var{entries} =} __xzip__ (@var{commandname}, @var{extension}, @var{commandtemplate}, @var{files}, @var{outdir}) +## Undocumented internal function. +## @end deftypefn + +## Compress the list of files and/or directories specified in @var{files} +## with the external compression command @var{commandname}. The template +## @var{commandtemplate} is used to actually start the command. Each file +## is compressed separately and a new file with the extension @var{extension} +## is created and placed into the directory @var{outdir}. The original files +## are not touched. Existing compressed files are silently overwritten. +## This is an internal function. Do not use directly. + +function entries = __xzip__ (commandname, extension, + commandtemplate, files, outdir) + + if (nargin == 4 || nargin == 5) + if (! ischar (extension) || length (extension) == 0) + error ("__xzip__: extension has to be a string with finite length"); + endif + + if (nargin == 5 && ! exist (outdir, "dir")) + error ("__xzip__: output directory does not exist"); + endif + + if (ischar (files)) + files = cellstr (files); + else + error ("__xzip__: expecting FILES to be a character array"); + endif + + if (nargin == 4) + outdir = tmpnam (); + mkdir (outdir); + endif + + cwd = pwd(); + unwind_protect + files = glob (files); + + ## Ignore any file with the compress extension + files (cellfun (@(x) length(x) > length(extension) + && strcmp (x((end - length(extension) + 1):end), extension), + files)) = []; + + copyfile (files, outdir); + + [d, f] = myfileparts(files); + + cd (outdir); + + cmd = sprintf (commandtemplate, sprintf (" %s", f{:})); + + [status, output] = system (cmd); + if (status == 0) + + if (nargin == 5) + compressed_files = cellfun( + @(x) fullfile (outdir, sprintf ("%s.%s", x, extension)), + f, "UniformOutput", false); + else + movefile (cellfun(@(x) sprintf ("%s.%s", x, extension), f, + "UniformOutput", false), cwd); + ## FIXME this does not work when you try to compress directories + + compressed_files = cellfun(@(x) sprintf ("%s.%s", x, extension), + files, "UniformOutput", false); + endif + + if (nargout > 0) + entries = compressed_files; + endif + else + error ("%s command failed with exit status = %d", + commandname, status); + endif + unwind_protect_cleanup + cd(cwd); + if (nargin == 4) + crr = confirm_recursive_rmdir (); + unwind_protect + confirm_recursive_rmdir (false); + rmdir (outdir, "s"); + unwind_protect_cleanup + confirm_recursive_rmdir (crr); + end_unwind_protect + endif + end_unwind_protect + else + print_usage (); + endif + +endfunction + +function [d, f] = myfileparts (files) + [d, f, ext] = cellfun (@(x) fileparts (x), files, "UniformOutput", false); + f = cellfun (@(x, y) sprintf ("%s%s", x, y), f, ext, + "UniformOutput", false); + idx = cellfun (@(x) isdir (x), files); + d(idx) = ""; + f(idx) = files(idx); +endfunction + +## FIXME -- reinstate these tests if we invent a way to test private +## functions directly. +## +## %!error <extension has to be a string with finite length> +## %! __xzip__("gzip", "", "gzip -r %s", "bla"); +## %!error <no files to move> +## %! __xzip__("gzip", ".gz", "gzip -r %s", tmpnam); +## %!error <command failed with exit status> +## %! # test __xzip__ with invalid compression command +## %! unwind_protect +## %! filename = tmpnam; +## %! dummy = 1; +## %! save(filename, "dummy"); +## %! dirname = tmpnam; +## %! mkdir(dirname); +## %! entry = __xzip__("gzip", ".gz", "xxxzipxxx -r %s 2>/dev/null", +## %! filename, dirname); +## %! unwind_protect_cleanup +## %! delete(filename); +## %! rmdir(dirname); +## %! end_unwind_protect
--- a/scripts/optimization/__dogleg__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,63 +0,0 @@ -## Copyright (C) 2008, 2009 Jaroslav Hajek -## -## 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} {@var{x}} = __dogleg__ (@var{r}, @var{b}, @var{x}, @var{d}, @var{delta}) -## Undocumented internal function. -## @end deftypefn - -## Solve the double dogleg trust-region least-squares problem: -## Minimize norm(r*x-b) subject to the constraint norm(d.*x) <= delta, -## x being a convex combination of the gauss-newton and scaled gradient. - -## TODO: error checks -## TODO: handle singularity, or leave it up to mldivide? - -function x = __dogleg__ (r, b, d, delta) - ## Get Gauss-Newton direction. - x = r \ b; - xn = norm (d .* x); - if (xn > delta) - ## GN is too big, get scaled gradient. - s = (r' * b) ./ d; - sn = norm (s); - if (sn > 0) - ## Normalize and rescale. - s = (s / sn) ./ d; - ## Get the line minimizer in s direction. - tn = norm (r*s); - snm = (sn / tn) / tn; - if (snm < delta) - ## Get the dogleg path minimizer. - bn = norm (b); - dxn = delta/xn; snmd = snm/delta; - t = (bn/sn) * (bn/xn) * snmd; - t -= dxn * snmd^2 - sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2)); - alpha = dxn*(1-snmd^2) / t; - else - alpha = 0; - endif - else - alpha = delta / xn; - snm = 0; - endif - ## Form the appropriate convex combination. - x = alpha * x + ((1-alpha) * min (snm, delta)) * s; - endif -endfunction -
--- a/scripts/optimization/__doglegm__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,63 +0,0 @@ -## Copyright (C) 2008, 2009 Jaroslav Hajek -## -## 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} {@var{x}} = __doglegm__ (@var{r}, @var{b}, @var{x}, @var{d}, @var{delta}) -## @end deftypefn - -## Solve the double dogleg trust-region minimization problem: -## Minimize 1/2*norm(r*x)^2 subject to the constraint norm(d.*x) <= delta, -## x being a convex combination of the gauss-newton and scaled gradient. - -## TODO: error checks -## TODO: handle singularity, or leave it up to mldivide? - -function x = __doglegm__ (r, g, d, delta) - ## Get Gauss-Newton direction. - b = r' \ g; - x = r \ b; - xn = norm (d .* x); - if (xn > delta) - ## GN is too big, get scaled gradient. - s = g ./ d; - sn = norm (s); - if (sn > 0) - ## Normalize and rescale. - s = (s / sn) ./ d; - ## Get the line minimizer in s direction. - tn = norm (r*s); - snm = (sn / tn) / tn; - if (snm < delta) - ## Get the dogleg path minimizer. - bn = norm (b); - dxn = delta/xn; snmd = snm/delta; - t = (bn/sn) * (bn/xn) * snmd; - t -= dxn * snmd^2 - sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2)); - alpha = dxn*(1-snmd^2) / t; - else - alpha = 0; - endif - else - alpha = delta / xn; - snm = 0; - endif - ## Form the appropriate convex combination. - x = alpha * x + ((1-alpha) * min (snm, delta)) * s; - endif -endfunction -
--- a/scripts/optimization/__fdjac__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,48 +0,0 @@ -## Copyright (C) 2008, 2009 Jaroslav Hajek -## -## 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} {} __fdjac__ (@var{fcn}, @var{x}, @var{fvec}, @var{err}) -## Undocumented internal function. -## @end deftypefn - -function fjac = __fdjac__ (fcn, x, fvec, cdif, err = 0) - if (cdif) - err = (max (eps, err)) ^ (1/3); - h = max (abs (x), 1)*err; # FIXME? - fjac = zeros (length (fvec), numel (x)); - for i = 1:numel (x) - x1 = x2 = x; - x1(i) += h(i); - x2(i) -= h(i); - fjac(:,i) = (fcn (x1)(:) - fcn (x2)(:)) / (x1(i) - x2(i)); - endfor - else - err = sqrt (max (eps, err)); - h = max (abs (x), 1)*err; # FIXME? - fjac = zeros (length (fvec), numel (x)); - for i = 1:numel (x) - x1 = x; - x1(i) += h(i); - fjac(:,i) = (fcn (x1)(:) - fvec) / (x1(i) - x(i)); - endfor - endif -endfunction - - -
--- a/scripts/optimization/fminunc.m Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/optimization/fminunc.m Tue Dec 01 22:40:37 2009 -0500 @@ -349,3 +349,43 @@ %! assert (x, ones (1, 4), tol); %! assert (fval, 0, tol); +## Solve the double dogleg trust-region minimization problem: +## Minimize 1/2*norm(r*x)^2 subject to the constraint norm(d.*x) <= delta, +## x being a convex combination of the gauss-newton and scaled gradient. + +## TODO: error checks +## TODO: handle singularity, or leave it up to mldivide? + +function x = __doglegm__ (r, g, d, delta) + ## Get Gauss-Newton direction. + b = r' \ g; + x = r \ b; + xn = norm (d .* x); + if (xn > delta) + ## GN is too big, get scaled gradient. + s = g ./ d; + sn = norm (s); + if (sn > 0) + ## Normalize and rescale. + s = (s / sn) ./ d; + ## Get the line minimizer in s direction. + tn = norm (r*s); + snm = (sn / tn) / tn; + if (snm < delta) + ## Get the dogleg path minimizer. + bn = norm (b); + dxn = delta/xn; snmd = snm/delta; + t = (bn/sn) * (bn/xn) * snmd; + t -= dxn * snmd^2 - sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2)); + alpha = dxn*(1-snmd^2) / t; + else + alpha = 0; + endif + else + alpha = delta / xn; + snm = 0; + endif + ## Form the appropriate convex combination. + x = alpha * x + ((1-alpha) * min (snm, delta)) * s; + endif +endfunction
--- a/scripts/optimization/fsolve.m Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/optimization/fsolve.m Tue Dec 01 22:40:37 2009 -0500 @@ -526,3 +526,43 @@ %! assert (norm (f) < tol); %! assert (norm (x - x_opt, Inf) < tol); +## Solve the double dogleg trust-region least-squares problem: +## Minimize norm(r*x-b) subject to the constraint norm(d.*x) <= delta, +## x being a convex combination of the gauss-newton and scaled gradient. + +## TODO: error checks +## TODO: handle singularity, or leave it up to mldivide? + +function x = __dogleg__ (r, b, d, delta) + ## Get Gauss-Newton direction. + x = r \ b; + xn = norm (d .* x); + if (xn > delta) + ## GN is too big, get scaled gradient. + s = (r' * b) ./ d; + sn = norm (s); + if (sn > 0) + ## Normalize and rescale. + s = (s / sn) ./ d; + ## Get the line minimizer in s direction. + tn = norm (r*s); + snm = (sn / tn) / tn; + if (snm < delta) + ## Get the dogleg path minimizer. + bn = norm (b); + dxn = delta/xn; snmd = snm/delta; + t = (bn/sn) * (bn/xn) * snmd; + t -= dxn * snmd^2 - sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2)); + alpha = dxn*(1-snmd^2) / t; + else + alpha = 0; + endif + else + alpha = delta / xn; + snm = 0; + endif + ## Form the appropriate convex combination. + x = alpha * x + ((1-alpha) * min (snm, delta)) * s; + endif +endfunction +
--- a/scripts/optimization/module.mk Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/optimization/module.mk Tue Dec 01 22:40:37 2009 -0500 @@ -1,21 +1,22 @@ FCN_FILE_DIRS += optimization +optimization_PRIVATE_FCN_FILES = \ + optimization/private/__fdjac__.m + optimization_FCN_FILES = \ + optimization/__all_opts__.m \ + optimization/fminunc.m \ + optimization/fsolve.m \ optimization/fzero.m \ - optimization/__fdjac__.m \ - optimization/__dogleg__.m \ - optimization/__doglegm__.m \ - optimization/fsolve.m \ - optimization/fminunc.m \ optimization/glpk.m \ optimization/glpkmex.m \ optimization/lsqnonneg.m \ - optimization/pqpnonneg.m \ + optimization/optimget.m \ optimization/optimset.m \ - optimization/optimget.m \ - optimization/__all_opts__.m \ + optimization/pqpnonneg.m \ optimization/qp.m \ - optimization/sqp.m + optimization/sqp.m \ + $(optimization_PRIVATE_FCN_FILES) FCN_FILES += $(optimization_FCN_FILES)
--- a/scripts/path/__extractpath__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,95 +0,0 @@ -## Copyright (C) 2005, 2006, 2007, 2008, 2009 Bill Denney -## Copyright (C) 2007 Ben Abbott -## -## 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} {@var{val} =} __extractpath__ (@var{file}) -## Undocumented internal function. -## @end deftypefn - -## Extact the path information from the script/function @var{file}, -## created by @file{savepath.m}. If @var{file} is omitted, -## @file{~/.octaverc} is used. If successful, @code{__extractpath__} -## returns the path specified in @var{file}. - -## Author: Ben Abbott <bpabbott@mac.com> - -function specifiedpath = __extractpath__ (savefile) - - ## The majority of this code was borrowed from savepath.m. - ## FIXME -- is there some way to share the common parts instead of - ## duplicating? - - beginstring = "## Begin savepath auto-created section, do not edit"; - endstring = "## End savepath auto-created section"; - - if (nargin == 0) - savefile = tilde_expand ("~/.octaverc"); - endif - - ## Parse the file if it exists to see if we should replace a section - ## or create a section. - startline = 0; - endline = 0; - filelines = {}; - if (exist (savefile) == 2) - ## read in all lines of the file - [fid, msg] = fopen (savefile, "rt"); - if (fid < 0) - error ("__extractpath__: could not open savefile, %s: %s", savefile, msg); - endif - unwind_protect - linenum = 0; - while (linenum >= 0) - result = fgetl (fid); - if (isnumeric (result)) - ## End at the end of file. - linenum = -1; - else - linenum++; - filelines{linenum} = result; - ## Find the first and last lines if they exist in the file. - if (strcmp (result, beginstring)) - startline = linenum + 1; - elseif (strcmp (result, endstring)) - endline = linenum - 1; - endif - endif - endwhile - unwind_protect_cleanup - closeread = fclose (fid); - if (closeread < 0) - error ("savepath: could not close savefile after reading, %s", - savefile); - endif - end_unwind_protect - endif - - ## Extract the path specifiation. - if (startline > endline || (startline > 0 && endline == 0)) - error ("savepath: unable to parse file, %s", savefile); - elseif (startline > 0) - ## Undo doubling of single quote characters performed by savepath. - specifiedpath = strrep (regexprep (cstrcat (filelines(startline:endline){:}), - " *path *\\('(.*)'\\); *", "$1"), - "''", "'"); - else - specifiedpath = ""; - endif - -endfunction
--- a/scripts/path/module.mk Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/path/module.mk Tue Dec 01 22:40:37 2009 -0500 @@ -1,7 +1,6 @@ FCN_FILE_DIRS += path path_FCN_FILES = \ - path/__extractpath__.m \ path/matlabroot.m \ path/pathdef.m \ path/savepath.m
--- a/scripts/path/pathdef.m Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/path/pathdef.m Tue Dec 01 22:40:37 2009 -0500 @@ -1,4 +1,5 @@ -## Copyright (C) 2008, 2009 Ben Abbott +## Copyright (C) 2005, 2006, 2007, 2008, 2009 Bill Denney +## Copyright (C) 2007, 2008, 2009 Ben Abbott ## ## This file is part of Octave. ## @@ -59,3 +60,75 @@ endif endfunction + +## Extact the path information from the script/function @var{file}, +## created by @file{savepath.m}. If @var{file} is omitted, +## @file{~/.octaverc} is used. If successful, @code{__extractpath__} +## returns the path specified in @var{file}. + +## Author: Ben Abbott <bpabbott@mac.com> + +function specifiedpath = __extractpath__ (savefile) + + ## The majority of this code was borrowed from savepath.m. + ## FIXME -- is there some way to share the common parts instead of + ## duplicating? + + beginstring = "## Begin savepath auto-created section, do not edit"; + endstring = "## End savepath auto-created section"; + + if (nargin == 0) + savefile = tilde_expand ("~/.octaverc"); + endif + + ## Parse the file if it exists to see if we should replace a section + ## or create a section. + startline = 0; + endline = 0; + filelines = {}; + if (exist (savefile) == 2) + ## read in all lines of the file + [fid, msg] = fopen (savefile, "rt"); + if (fid < 0) + error ("__extractpath__: could not open savefile, %s: %s", savefile, msg); + endif + unwind_protect + linenum = 0; + while (linenum >= 0) + result = fgetl (fid); + if (isnumeric (result)) + ## End at the end of file. + linenum = -1; + else + linenum++; + filelines{linenum} = result; + ## Find the first and last lines if they exist in the file. + if (strcmp (result, beginstring)) + startline = linenum + 1; + elseif (strcmp (result, endstring)) + endline = linenum - 1; + endif + endif + endwhile + unwind_protect_cleanup + closeread = fclose (fid); + if (closeread < 0) + error ("savepath: could not close savefile after reading, %s", + savefile); + endif + end_unwind_protect + endif + + ## Extract the path specifiation. + if (startline > endline || (startline > 0 && endline == 0)) + error ("savepath: unable to parse file, %s", savefile); + elseif (startline > 0) + ## Undo doubling of single quote characters performed by savepath. + specifiedpath = strrep (regexprep (cstrcat (filelines(startline:endline){:}), + " *path *\\('(.*)'\\); *", "$1"), + "''", "'"); + else + specifiedpath = ""; + endif + +endfunction
--- a/scripts/plot/__add_datasource__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,55 +0,0 @@ -## Copyright (C) 2008, 2009 David Bateman -## -## 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} {@var{newargs} =} __add_datasource__ (@var{fcn}, @var{h}, @var{data}, @var{varargin}) -## Undocumented internal function. -## @end deftypefn - -function newargs = __add_datasource__ (fcn, h, data, varargin) - - if (nargin < 3) - error ("internal error"); - endif - - if (ischar (data)) - data = {data}; - endif - - for i = 1 : numel (data) - addproperty (strcat (data{i}, "datasource"), h, "string", ""); - endfor - - i = 0; - newargs = {}; - while (i < numel (varargin)) - arg = varargin{++i}; - if (i != numel(varargin) && ischar (arg) - && length (arg) > 9 && strcmpi (arg(end-9:end), "datasource")) - arg = tolower (arg); - val = varargin{++i}; - if (ischar (val)) - set (h, arg, val); - else - error ("%s: expecting data source to be a string", fcn); - endif - else - newargs{end + 1} = arg; - endif - endwhile -endfunction
--- a/scripts/statistics/base/__quantile__.m Tue Dec 01 16:21:12 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,262 +0,0 @@ -## Copyright (C) 2008, 2009 Ben Abbott and Jaroslav Hajek -## -## 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} {@var{q} =} __quantile__ (@var{x}, @var{p}) -## @deftypefnx {Function File} {@var{q} =} __quantile__ (@var{x}, @var{p}, @var{method}) -## Undocumented internal function. -## @end deftypefn - -## For the cumulative probability values in @var{p}, compute the -## quantiles, @var{q} (the inverse of the cdf), for the sample, @var{x}. -## -## The optional input, @var{method}, refers to nine methods available in R -## (http://www.r-project.org/). The default is @var{method} = 7. For more -## detail, see `help quantile'. -## @seealso{prctile, quantile, statistics} - -## Author: Ben Abbott <bpabbott@mac.com> -## Vectorized version: Jaroslav Hajek <highegg@gmail.com> -## Description: Quantile function of a empirical samples - -function inv = __quantile__ (x, p, method = 5) - - if (nargin < 2 || nargin > 3) - print_usage (); - endif - - if (! ismatrix (x)) - error ("quantile: x must be a matrix"); - endif - - ## Save length and set shape of quantiles. - n = numel (p); - p = p(:); - - ## Save length and set shape of samples. - ## FIXME: does sort guarantee that NaN's come at the end? - x = sort (x); - m = sum (! isnan (x)); - mx = size (x, 1); - nx = size (x, 2); - - ## Initialize output values. - inv = Inf*(-(p < 0) + (p > 1)); - inv = repmat (inv, 1, nx); - - ## Do the work. - if (any(k = find((p >= 0) & (p <= 1)))) - n = length (k); - p = p (k); - ## Special case. - if (mx == 1) - inv(k,:) = repmat (x, n, 1); - return - endif - - ## The column-distribution indices. - pcd = kron (ones (n, 1), mx*(0:nx-1)); - mm = kron (ones (n, 1), m); - switch method - case {1, 2, 3} - switch method - case 1 - p = max (ceil (kron (p, m)), 1); - inv(k,:) = x(p + pcd); - - case 2 - p = kron (p, m); - p_lr = max (ceil (p), 1); - p_rl = min (floor (p + 1), mm); - inv(k,:) = (x(p_lr + pcd) + x(p_rl + pcd))/2; - - case 3 - ## Used by SAS, method PCTLDEF=2. - ## http://support.sas.com/onlinedoc/913/getDoc/en/statug.hlp/stdize_sect14.htm - t = max (kron (p, m), 1); - t = roundb (t); - inv(k,:) = x(t + pcd); - endswitch - - otherwise - switch method - case 4 - p = kron (p, m); - - case 5 - ## Used by Matlab. - p = kron (p, m) + 0.5; - - case 6 - ## Used by Minitab and SPSS. - p = kron (p, m+1); - - case 7 - ## Used by S and R. - p = kron (p, m-1) + 1; - - case 8 - ## Median unbiased . - p = kron (p, m+1/3) + 1/3; - - case 9 - ## Approximately unbiased respecting order statistics. - p = kron (p, m+0.25) + 0.375; - - otherwise - error ("quantile: Unknown method, '%d'", method); - endswitch - - ## Duplicate single values. - imm1 = mm == 1; - x(2,imm1) = x(1,imm1); - - ## Interval indices. - pi = max (min (floor (p), mm-1), 1); - pr = max (min (p - pi, 1), 0); - pi += pcd; - inv(k,:) = (1-pr) .* x(pi) + pr .* x(pi+1); - endswitch - endif - -endfunction - -%!test -%! p = 0.5; -%! x = sort (rand (11)); -%! q = __quantile__ (x, p); -%! assert (q, x(6,:)) - -%!test -%! p = [0.00, 0.25, 0.50, 0.75, 1.00]; -%! x = [1; 2; 3; 4]; -%! a = [1.0000 1.0000 2.0000 3.0000 4.0000 -%! 1.0000 1.5000 2.5000 3.5000 4.0000 -%! 1.0000 1.0000 2.0000 3.0000 4.0000 -%! 1.0000 1.0000 2.0000 3.0000 4.0000 -%! 1.0000 1.5000 2.5000 3.5000 4.0000 -%! 1.0000 1.2500 2.5000 3.7500 4.0000 -%! 1.0000 1.7500 2.5000 3.2500 4.0000 -%! 1.0000 1.4167 2.5000 3.5833 4.0000 -%! 1.0000 1.4375 2.5000 3.5625 4.0000]; -%! for m = (1:9) -%! q = __quantile__ (x, p, m).'; -%! assert (q, a(m,:), 0.0001) -%! endfor - -%!test -%! p = [0.00, 0.25, 0.50, 0.75, 1.00]; -%! x = [1; 2; 3; 4; 5]; -%! a = [1.0000 2.0000 3.0000 4.0000 5.0000 -%! 1.0000 2.0000 3.0000 4.0000 5.0000 -%! 1.0000 1.0000 2.0000 4.0000 5.0000 -%! 1.0000 1.2500 2.5000 3.7500 5.0000 -%! 1.0000 1.7500 3.0000 4.2500 5.0000 -%! 1.0000 1.5000 3.0000 4.5000 5.0000 -%! 1.0000 2.0000 3.0000 4.0000 5.0000 -%! 1.0000 1.6667 3.0000 4.3333 5.0000 -%! 1.0000 1.6875 3.0000 4.3125 5.0000]; -%! for m = (1:9) -%! q = __quantile__ (x, p, m).'; -%! assert (q, a(m,:), 0.0001) -%! endfor - -%!test -%! p = [0.00, 0.25, 0.50, 0.75, 1.00]; -%! x = [1; 2; 5; 9]; -%! a = [1.0000 1.0000 2.0000 5.0000 9.0000 -%! 1.0000 1.5000 3.5000 7.0000 9.0000 -%! 1.0000 1.0000 2.0000 5.0000 9.0000 -%! 1.0000 1.0000 2.0000 5.0000 9.0000 -%! 1.0000 1.5000 3.5000 7.0000 9.0000 -%! 1.0000 1.2500 3.5000 8.0000 9.0000 -%! 1.0000 1.7500 3.5000 6.0000 9.0000 -%! 1.0000 1.4167 3.5000 7.3333 9.0000 -%! 1.0000 1.4375 3.5000 7.2500 9.0000]; -%! for m = (1:9) -%! q = __quantile__ (x, p, m).'; -%! assert (q, a(m,:), 0.0001) -%! endfor - -%!test -%! p = [0.00, 0.25, 0.50, 0.75, 1.00]; -%! x = [1; 2; 5; 9; 11]; -%! a = [1.0000 2.0000 5.0000 9.0000 11.0000 -%! 1.0000 2.0000 5.0000 9.0000 11.0000 -%! 1.0000 1.0000 2.0000 9.0000 11.0000 -%! 1.0000 1.2500 3.5000 8.0000 11.0000 -%! 1.0000 1.7500 5.0000 9.5000 11.0000 -%! 1.0000 1.5000 5.0000 10.0000 11.0000 -%! 1.0000 2.0000 5.0000 9.0000 11.0000 -%! 1.0000 1.6667 5.0000 9.6667 11.0000 -%! 1.0000 1.6875 5.0000 9.6250 11.0000]; -%! for m = (1:9) -%! q = __quantile__ (x, p, m).'; -%! assert (q, a(m,:), 0.0001) -%! endfor - -%!test -%! p = [0.00, 0.25, 0.50, 0.75, 1.00]; -%! x = [16; 11; 15; 12; 15; 8; 11; 12; 6; 10]; -%! a = [6.0000 10.0000 11.0000 15.0000 16.0000 -%! 6.0000 10.0000 11.5000 15.0000 16.0000 -%! 6.0000 8.0000 11.0000 15.0000 16.0000 -%! 6.0000 9.0000 11.0000 13.5000 16.0000 -%! 6.0000 10.0000 11.5000 15.0000 16.0000 -%! 6.0000 9.5000 11.5000 15.0000 16.0000 -%! 6.0000 10.2500 11.5000 14.2500 16.0000 -%! 6.0000 9.8333 11.5000 15.0000 16.0000 -%! 6.0000 9.8750 11.5000 15.0000 16.0000]; -%! for m = (1:9) -%! q = __quantile__ (x, p, m).'; -%! assert (q, a(m,:), 0.0001) -%! endfor - -%!test -%! p = [0.00, 0.25, 0.50, 0.75, 1.00]; -%! x = [-0.58851; 0.40048; 0.49527; -2.551500; -0.52057; ... -%! -0.17841; 0.057322; -0.62523; 0.042906; 0.12337]; -%! a = [-2.551474 -0.588505 -0.178409 0.123366 0.495271 -%! -2.551474 -0.588505 -0.067751 0.123366 0.495271 -%! -2.551474 -0.625231 -0.178409 0.123366 0.495271 -%! -2.551474 -0.606868 -0.178409 0.090344 0.495271 -%! -2.551474 -0.588505 -0.067751 0.123366 0.495271 -%! -2.551474 -0.597687 -0.067751 0.192645 0.495271 -%! -2.551474 -0.571522 -0.067751 0.106855 0.495271 -%! -2.551474 -0.591566 -0.067751 0.146459 0.495271 -%! -2.551474 -0.590801 -0.067751 0.140686 0.495271]; -%! for m = (1:9) -%! q = __quantile__ (x, p, m).'; -%! assert (q, a(m,:), 0.0001) -%! endfor - -%!test -%! p = 0.5; -%! x = [0.112600, 0.114800, 0.052100, 0.236400, 0.139300 -%! 0.171800, 0.727300, 0.204100, 0.453100, 0.158500 -%! 0.279500, 0.797800, 0.329600, 0.556700, 0.730700 -%! 0.428800, 0.875300, 0.647700, 0.628700, 0.816500 -%! 0.933100, 0.931200, 0.963500, 0.779600, 0.846100]; -%! tol = 0.00001; -%! x(5,5) = NaN; -%! assert (__quantile__ (x, p), [0.27950, 0.79780, 0.32960, 0.55670, 0.44460], tol); -%! x(1,1) = NaN; -%! assert (__quantile__ (x, p), [0.35415, 0.79780, 0.32960, 0.55670, 0.44460], tol); -%! x(3,3) = NaN; -%! assert (__quantile__ (x, p), [0.35415, 0.79780, 0.42590, 0.55670, 0.44460], tol); -
--- a/scripts/statistics/base/module.mk Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/statistics/base/module.mk Tue Dec 01 22:40:37 2009 -0500 @@ -1,7 +1,6 @@ FCN_FILE_DIRS += statistics/base statistics_base_FCN_FILES = \ - statistics/base/__quantile__.m \ statistics/base/center.m \ statistics/base/cloglog.m \ statistics/base/cor.m \
--- a/scripts/statistics/base/quantile.m Tue Dec 01 16:21:12 2009 -0500 +++ b/scripts/statistics/base/quantile.m Tue Dec 01 22:40:37 2009 -0500 @@ -1,4 +1,4 @@ -## Copyright (C) 2008, 2009 Ben Abbott +## Copyright (C) 2008, 2009 Ben Abbott and Jaroslav Hajek ## ## This file is part of Octave. ## @@ -265,3 +265,116 @@ %! yexp = median (x, dim); %! assert (yobs, yexp); +## For the cumulative probability values in @var{p}, compute the +## quantiles, @var{q} (the inverse of the cdf), for the sample, @var{x}. +## +## The optional input, @var{method}, refers to nine methods available in R +## (http://www.r-project.org/). The default is @var{method} = 7. For more +## detail, see `help quantile'. +## @seealso{prctile, quantile, statistics} + +## Author: Ben Abbott <bpabbott@mac.com> +## Vectorized version: Jaroslav Hajek <highegg@gmail.com> +## Description: Quantile function of a empirical samples + +function inv = __quantile__ (x, p, method = 5) + + if (nargin < 2 || nargin > 3) + print_usage (); + endif + + if (! ismatrix (x)) + error ("quantile: x must be a matrix"); + endif + + ## Save length and set shape of quantiles. + n = numel (p); + p = p(:); + + ## Save length and set shape of samples. + ## FIXME: does sort guarantee that NaN's come at the end? + x = sort (x); + m = sum (! isnan (x)); + mx = size (x, 1); + nx = size (x, 2); + + ## Initialize output values. + inv = Inf*(-(p < 0) + (p > 1)); + inv = repmat (inv, 1, nx); + + ## Do the work. + if (any(k = find((p >= 0) & (p <= 1)))) + n = length (k); + p = p (k); + ## Special case. + if (mx == 1) + inv(k,:) = repmat (x, n, 1); + return + endif + + ## The column-distribution indices. + pcd = kron (ones (n, 1), mx*(0:nx-1)); + mm = kron (ones (n, 1), m); + switch method + case {1, 2, 3} + switch method + case 1 + p = max (ceil (kron (p, m)), 1); + inv(k,:) = x(p + pcd); + + case 2 + p = kron (p, m); + p_lr = max (ceil (p), 1); + p_rl = min (floor (p + 1), mm); + inv(k,:) = (x(p_lr + pcd) + x(p_rl + pcd))/2; + + case 3 + ## Used by SAS, method PCTLDEF=2. + ## http://support.sas.com/onlinedoc/913/getDoc/en/statug.hlp/stdize_sect14.htm + t = max (kron (p, m), 1); + t = roundb (t); + inv(k,:) = x(t + pcd); + endswitch + + otherwise + switch method + case 4 + p = kron (p, m); + + case 5 + ## Used by Matlab. + p = kron (p, m) + 0.5; + + case 6 + ## Used by Minitab and SPSS. + p = kron (p, m+1); + + case 7 + ## Used by S and R. + p = kron (p, m-1) + 1; + + case 8 + ## Median unbiased . + p = kron (p, m+1/3) + 1/3; + + case 9 + ## Approximately unbiased respecting order statistics. + p = kron (p, m+0.25) + 0.375; + + otherwise + error ("quantile: Unknown method, '%d'", method); + endswitch + + ## Duplicate single values. + imm1 = mm == 1; + x(2,imm1) = x(1,imm1); + + ## Interval indices. + pi = max (min (floor (p), mm-1), 1); + pr = max (min (p - pi, 1), 0); + pi += pcd; + inv(k,:) = (1-pr) .* x(pi) + pr .* x(pi+1); + endswitch + endif + +endfunction