Mercurial > hg > octave-max
changeset 9899:9f25290a35e8
more private function and subfunction changes
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -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 +++ b/scripts/general/isequal.m @@ -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 +++ b/scripts/general/isequalwithequalnans.m @@ -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 +++ b/scripts/general/module.mk @@ -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)
rename from scripts/general/__isequal__.m rename to scripts/general/private/__isequal__.m --- a/scripts/general/__isequal__.m +++ b/scripts/general/private/__isequal__.m @@ -181,50 +181,3 @@ 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/help/module.mk +++ b/scripts/help/module.mk @@ -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)
rename from scripts/help/__additional_help_message__.m rename to scripts/help/private/__additional_help_message__.m
deleted file mode 100644 --- a/scripts/image/__img__.m +++ /dev/null @@ -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
deleted file mode 100644 --- a/scripts/image/__img_via_file__.m +++ /dev/null @@ -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 +++ b/scripts/image/image.m @@ -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 +++ b/scripts/image/image_viewer.m @@ -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 +++ b/scripts/image/module.mk @@ -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/module.mk +++ b/scripts/miscellaneous/module.mk @@ -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)
rename from scripts/miscellaneous/__xzip__.m rename to scripts/miscellaneous/private/__xzip__.m --- a/scripts/miscellaneous/__xzip__.m +++ b/scripts/miscellaneous/private/__xzip__.m @@ -120,21 +120,24 @@ 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 +## 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
deleted file mode 100644 --- a/scripts/optimization/__dogleg__.m +++ /dev/null @@ -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 -
deleted file mode 100644 --- a/scripts/optimization/__doglegm__.m +++ /dev/null @@ -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 -
deleted file mode 100644 --- a/scripts/optimization/__fdjac__.m +++ /dev/null @@ -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 +++ b/scripts/optimization/fminunc.m @@ -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 +++ b/scripts/optimization/fsolve.m @@ -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 +++ b/scripts/optimization/module.mk @@ -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)
deleted file mode 100644 --- a/scripts/path/__extractpath__.m +++ /dev/null @@ -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 +++ b/scripts/path/module.mk @@ -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 +++ b/scripts/path/pathdef.m @@ -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
deleted file mode 100644 --- a/scripts/plot/__add_datasource__.m +++ /dev/null @@ -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
deleted file mode 100644 --- a/scripts/statistics/base/__quantile__.m +++ /dev/null @@ -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 +++ b/scripts/statistics/base/module.mk @@ -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 +++ b/scripts/statistics/base/quantile.m @@ -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