# HG changeset patch
# User Nicholas R. Jankowski <jankowskin@asme.org>
# Date 1607715939 18000
#      Fri Dec 11 14:45:39 2020 -0500
# Node ID e0d43c3f9c5284477a232efdb542ea8a08f1d793
# Parent  7d262c7f8b1d30457d06f021245356df3054fa8e
iqr quantile method change, extend dimension inputs, and speed increase (bug #59636)

* iqr.m: quantile method now matlab compatible, add "all" and vector dimensions. Update tests.
* NEWS: Note change.

diff -r 7d262c7f8b1d -r e0d43c3f9c52 NEWS
--- a/NEWS	Tue Dec 08 23:12:15 2020 -0500
+++ b/NEWS	Fri Dec 11 14:45:39 2020 -0500
@@ -179,6 +179,11 @@
 - The functions `quantile` and `prctile` now permit operating on dimensions
 greater than `ndims(x)`.
 
+- The function `iqr` now uses Matlab compatible interpolation for quantile
+values, demension input has been extended to include vector dimensions,
+"all", and dimensions greater than `ndims(x)`, and compatibly handles `Inf` and
+`NaN` input values.
+
 ### Alphabetical list of new functions added in Octave 7
 
 * `cospi`
diff -r 7d262c7f8b1d -r e0d43c3f9c52 scripts/statistics/iqr.m
--- a/scripts/statistics/iqr.m	Tue Dec 08 23:12:15 2020 -0500
+++ b/scripts/statistics/iqr.m	Fri Dec 11 14:45:39 2020 -0500
@@ -1,103 +1,278 @@
-########################################################################
-##
-## Copyright (C) 1995-2020 The Octave Project Developers
-##
-## See the file COPYRIGHT.md in the top-level directory of this
-## distribution or <https://octave.org/copyright/>.
-##
-## 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
-## <https://www.gnu.org/licenses/>.
-##
-########################################################################
-
-## -*- texinfo -*-
-## @deftypefn  {} {} iqr (@var{x})
-## @deftypefnx {} {} iqr (@var{x}, @var{dim})
-## Return the interquartile range, i.e., the difference between the upper
-## and lower quartile of the input data.
-##
-## If @var{x} is a matrix, do the above for first non-singleton dimension of
-## @var{x}.
-##
-## If the optional argument @var{dim} is given, operate along this dimension.
-##
-## As a measure of dispersion, the interquartile range is less affected by
-## outliers than either @code{range} or @code{std}.
-## @seealso{bounds, mad, range, std}
-## @end deftypefn
-
-function y = iqr (x, dim)
-
-  if (nargin < 1)
+########################################################################
+##
+## Copyright (C) 1996-2020 The Octave Project Developers
+##
+## See the file COPYRIGHT.md in the top-level directory of this
+## distribution or <https://octave.org/copyright/>.
+##
+## 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
+## <https://www.gnu.org/licenses/>.
+##
+########################################################################
+
+## -*- texinfo -*-
+## @deftypefn {} {@var{Z} =} iqr (@var{X})
+## @deftypefnx {} {@var{Z} =} iqr (@var{X}, @var{DIM})
+## @deftypefnx {} {@var{Z} =} iqr (@var{X}, @qcode{"ALL"})
+## Return the interquartile range of @var{X}, defined as the distance between
+## the 25th and 75th percetile values of @var{X} calculated using:
+##    quantile (x, [0.25 0.75])
+##
+## If @var{X} is a vector, @code{iqr (@var{X})} will operate on the data in
+## @var{X}.
+##
+## If @var{X} is a matrix, @code{iqr (@var{X})} will operate independently on
+## each column in @var{X} returning a row vector @var{Z}.
+##
+## If @var{X} is a n-dimensional array, @code{iqr (@var{X})} will operate
+## independently on the first non-singleton dimension in @var{X}, returning an
+## array @var{Z} the same shape as @var{X} with the non-singleton dimenion
+## reduced to 1.
+##
+## The optional variable @var{DIM} can be used to force @code{iqr} to operate
+## over the specified dimension.  @var{DIM} can either be a scalar dimension or
+## a vector of non-repeating dimensions over which to operater.  In either case
+## @var{DIM} must be positive integers. A vector @var{DIM} concatenates all
+## specified deminsions for independent operation by @code{iqr}.
+##
+## Specifying dimension @qcode{"ALL"} will force @code{iqr} to operate
+## on all elements of @var{X}, and is equivalent to @code{iqr (@var{X}(:))}.
+## Similarly, specifying a vector dimension including all non-singleton
+## dimensions of @var{X} is equivalent to @code{iqr (@var{X}, @qcode{"ALL"})}.
+##
+## If @var{X} is a scalar, or only singleton dimensions are specified for
+## @var{DIM}, the output will be @code{zeros (size (@var{X}))}.
+##
+## @seealso{mad, range, prctile, quantile}
+## @end deftypefn
+
+## Author: Nicholas Jankowski <jankowskin@asme.org>
+## Created: 2020-11-25
+
+function z = iqr (varargin)
+  ##TODO:  when octave develops Probability Distribution Obects, enable
+  ## Matlab compatible iqr handling for those object types.
+
+  ## input checks
+  if (nargin == 1)
+    x = varargin{1};
+    dim = [];
+  elseif (nargin == 2)
+    x = varargin{1};
+    dim = varargin{2};
+  else
     print_usage ();
   endif
 
-  if (! (isnumeric (x) || islogical (x)))
-    error ("iqr: X must be a numeric vector or matrix");
+  if (~ (isnumeric (x) || islogical (x)))
+    error ("iqr: X must be a numeric vector or matrix");
   endif
+
+  vecdim_flag = false;
+  sz_x = size (x);
+  x_dims = ndims (x);
+
+  if isempty (dim)
+    ## find first non-singleton dim, works for n-dimensional vectors, too.
+    if (max (sz_x) == 1)
+      dim = 2;
+    else
+      dim = find ((sz_x > 1), 1);
+    endif
+  else
+
+    if (isnumeric (dim) && all (dim > 0) && all (rem (dim,1) == 0))
+      if (isscalar (dim))
+        ##nothing to do, skip other checks/errors.
+      elseif (isvector (dim) && ((num_vecdims = numel (dim)) > 1) ...
+        && all (diff (sort (dim))))
+        ## vecdim must be 1D, non repeating
+
+        ##catch trivial case of vecdim being all dimensions (same as "all")
+        highest_dim = (max (x_dims, max (dim)));
+        if (num_vecdims == x_dims) && (highest_dim == x_dims)
+          x = x(:);
+          sz_x = size (x);
+          dim = 1;
+
+        else
+        ## move dimensions for operation to the front, keeping the rest in order
+        ## reshape those dims into dim1 (column) vectors
+        ## process as normal for a dim1 iqr on X, reshape when done.
+
+          vecdim_flag = true;  ## flag for final reshape
+
+          ## allow row or column dims, process as row vector
+          if iscolumn (dim)
+            dim = dim';
+          endif
+
+          ##set the permutation vector with dim at front and rest trailing,
+          perm = [1 : highest_dim];
+          perm(dim)=[];
+          perm = [dim, perm];
+
+          ##reshape x to put dims to process at front
+          x = permute(x, perm);
+          sz_x_new = size (x);
+
+          ## preserve trailing singletons when dim > ndims (x)
+          sz_x_new = [sz_x_new, ones(1, highest_dim - numel (sz_x_new))];
+
+          newshape = [prod(sz_x_new(1:num_vecdims)), ones(1,(num_vecdims-1)), sz_x_new((num_vecdims+1):end)];
+
+          if numel(newshape) == 1
+            newshape = [newshape, 1];
+          endif
+
+          ##collapse process dimensions into single column for iqr calc
+          x = reshape (x, newshape);
+
+          ##reset dim to operate on new column vectors
+          dim = 1;
+        endif
+
+      else
+        error ("iqr: vector dimension must be non-repeating positive integer vector");
+      endif
+
+    elseif (strcmp (tolower (dim), "all"))
+      ## "ALL" simplifies to collapsing all elements to single vector
+      x = x(:);
+      dim = 1;
+      sz_x = size (x);
+
+    else
+      error ("iqr: dimension must be a positive integer scalar, vector, or 'all'");
+    endif
+
+  endif
+
+  if (((dim > x_dims) || (sz_x(dim) == 1)) && all (isfinite(x)))
+    ## shortcut easy zeros
+    z = zeros (sz_x);
+
+  elseif (iscolumn (x) && (dim == 1))
+    ## catch col vector with quantile/diff dim requirement mismatch
+    z = abs (diff (quantile (x, [0.25 0.75], 1), [], 2));
+
+  else
+    z = abs (diff (quantile (x, [0.25 0.75], dim), [], dim));
+  endif
+
+  if (vecdim_flag)
+    z = ipermute(z, perm);
+  endif
+
+endfunction
 
-  nd = ndims (x);
-  sz = size (x);
-  nel = numel (x);
-  if (nargin != 2)
-    ## Find the first non-singleton dimension.
-    (dim = find (sz > 1, 1)) || (dim = 1);
-  else
-    if (!(isscalar (dim) && dim == fix (dim))
-        || !(1 <= dim && dim <= nd))
-      error ("iqr: DIM must be an integer and a valid dimension");
-    endif
-  endif
+%!assert (iqr (17), 0)
+%!assert (iqr (17, 1), 0)
+%!assert (iqr (17, 4), 0)
+%!assert (iqr ([1:3]), 1.5)
+%!assert (iqr ([1:4]), 2)
+%!assert (iqr ([1:5]), 2.5)
+%!assert (iqr ([1:10]), 5)
+%!assert (iqr ([1:10]'), 5)
+%!assert (iqr ([1:10],2), 5)
+%!assert (iqr ([1:10],1), zeros(1,10))
+%!assert (iqr ([1:10],3), zeros(1,10))
+%!assert (iqr ([[1:5];[2:6]], "all"), 3)
+
+%!test
+%! x = reshape([1:6], [1 2 3]);
+%! assert (iqr(x), ones(1,1,3));
+%! assert (iqr(x,1), zeros (1,2,3));
+%! assert (iqr(x,2), ones(1,1,3));
+%! assert (iqr(x,3), [3 3]);
+
+## nD array tests
+%!test
+%! x = magic (4); x = cat (3,x, 2*x, 3*x); x = cat (4, x, 2*x);
+%! y = cat (3, 8*[1 1 1 1], 16*[1 1 1 1], 24*[1 1 1 1]);
+%! assert (iqr (x), cat (4, y, 2*y));
+%! assert (iqr (x, 1), cat (4, y, 2*y));
+%! y = cat (3, 4*[3 1 1 3]', 8*[3 1 1 3]', 12*[3 1 1 3]');
+%! assert (iqr (x, 2), cat (4, y, 2*y));
+%! y = [24 3 4.5 19.5; 7.5 16.5 15 12; 13.5 10.5 9, 18; 6 21 22.5 1.5];
+%! assert (iqr (x, 3), cat (4, y, 2*y));
+%! y = [16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
+%! assert (iqr (x, 4), cat (3, y, 2*y, 3*y));
+%! assert (iqr (x, 5), zeros (size (x)));
+
+## Vector dimension tests
+%!assert (iqr (17, [1 8]), 0)
+%!assert (iqr ([[1 2 5]; [2 5 6]], [1 2]), 3)
+%!assert (iqr (cat (3, [1 2 5; 2 5 6], [1 2 5; 2 5 6]), [1 2]), cat(3, 3, 3))
+%!assert (iqr (cat (3, [1 2 5; 2 5 6], [1 2 5; 2 5 6]), [1 2]'), cat(3, 3, 3))
+%!test
+%! x = magic (4); x = cat (3,x, 2*x, 3*x); x = cat (4, x, 2*x);
+%! y = cat (3, 8, 16, 24);
+%! assert (iqr (x, [1 2]), cat (4, y, 2*y));
+%! y = [14, 18.5, 17.5 19.5];
+%! assert (iqr (x, [1 3]), cat (4, y, 2*y));
+%! y = [10.5 12.5 11.5 15.0000];
+%! assert (iqr (x, [1 4]), cat (3, y, 2*y, 3*y));
+%! assert (iqr (x, [1 5]), iqr (x, 1));
+%! y = [24 13 12 25.5]';
+%! assert (iqr (x, [2 3]), cat (4, y, 2*y));
+%! y = [17.5, 9, 8, 18.5]';
+%! assert (iqr (x, [2 4]), cat (3, y, 2*y, 3*y));
+%! assert (iqr (x, [3 4]), [32 4 6 26; 10 22 20 16; 18 14 12 24; 8 28 30 2]);
+%! assert (iqr (x, [3 4]), iqr (x, [4 3]));
+%! assert (iqr (x, [1 2 3]), cat (4, 17.5, 35));
+%! assert (iqr (x, [2 3 4]), [29.5 19.5 23 31]');
+%! assert (iqr (x, [1 3 4]), [22 28 22 30.5]);
+%! assert (iqr (x, [1 2 4]), cat (3, 11, 22, 33));
+%! assert (iqr (x, [1 2 5]), iqr (x, [1 2]));
+%! assert (iqr (x, [5 6]), zeros (size (x)));
+
+##Inf, NaN checks
+%!assert (iqr (Inf), NaN)
+%!assert (iqr (-Inf), NaN)
+%!assert (iqr (NaN), NaN)
+%!assert (iqr (NaN), NaN)
+%!assert (iqr ([1 2 Inf], 1), [0 0 NaN])
+%!assert (iqr ([1 2 Inf], 2), Inf)
+%!assert (iqr ([1 2 -Inf], 1), [0 0 NaN])
+%!assert (iqr ([1 2 -Inf], 2), Inf)
+%!assert (iqr ([1 2 3 NaN], 1), [0 0 0 NaN])
+%!assert (iqr ([1 2 3 NaN], 2), 1.5)
+%!assert (iqr ([1 NaN 2 3], 2), 1.5)
+%!assert (iqr (NaN(2), 1), [NaN NaN])
+%!assert (iqr (NaN(2), 2), [NaN NaN]')
+%!assert (iqr (NaN(2), 3),  NaN(2))
+%!assert (iqr ([[1 2 5], [2 NaN 6]], "all"), 3.5)
 
-  ## This code is a bit heavy, but is needed until empirical_inv
-  ## can take a matrix, rather than just a vector argument.
-  n = sz(dim);
-  sz(dim) = 1;
-  if (isa (x, "single"))
-    y = zeros (sz, "single");
-  else
-    y = zeros (sz);
-  endif
-  stride = prod (sz(1:dim-1));
-  for i = 1 : nel / n
-    offset = i;
-    offset2 = 0;
-    while (offset > stride)
-      offset -= stride;
-      offset2 += 1;
-    endwhile
-    offset += offset2 * stride * n;
-    rng = [0 : n-1] * stride + offset;
-
-    y(i) = diff (empirical_inv ([1/4, 3/4], x(rng)));
-  endfor
-
-endfunction
-
-
-%!assert (iqr (1:101), 50)
-%!assert (iqr (single (1:101)), single (50))
-
-## FIXME: iqr throws horrible error when running across a dimension that is 1.
-%!test
-%! x = [1:100]';
-%! assert (iqr (x, 1), 50);
-%! assert (iqr (x', 2), 50);
-
-%!error <Invalid call> iqr ()
-%!error iqr (1)
-%!error iqr (['A'; 'B'])
-%!error iqr (1:10, 3)
+## Test input validation
+%!error iqr ()
+%!error iqr (1, 2, 3)
+%!error <X must be a numeric> iqr (['A'; 'B'])
+%!error <dimension must be a positive integer> iqr (1, 'A')
+%!error <dimension must be a positive integer> iqr (1, 0)
+%!error <dimension must be a positive integer> iqr (1, -2)
+%!error <dimension must be a positive integer> iqr (1, 1.4)
+%!error <dimension must be a positive integer> iqr (1, [1 -2])
+%!error <dimension must be a positive integer> iqr (1, [1 1.4])
+%!error <dimension must be a positive integer> iqr ([1 2 3], NaN)
+%!error <dimension must be a positive integer> iqr ([1 2 3], [2 NaN])
+%!error <dimension must be a positive integer> iqr ([1 2 3], Inf)
+%!error <dimension must be a positive integer> iqr ([1 2 3], [2 Inf])
+%!error <vector dimension must be non-repeating> iqr ([1 2 3], [1 2 1])
+%!error <vector dimension must be .* vector> iqr (1, [1 2; 3 4])
+
+
+
