# HG changeset patch
# User Nicholas R. Jankowski <jankowskin@asme.org>
# Date 1507441373 14400
#      Sun Oct 08 01:42:53 2017 -0400
# Node ID 058f2dc7226e27f054958b1e8d0007d4c0bd7be8
# Parent  f80dc6db9d185d88f10171becb382ac74694d02e
integral2.m: update to add functional limits and tiled implementation (bug #52074)

* scripts/general/integral2.m: added both functional y integration limits and david's 'tiled' method.

diff -r f80dc6db9d18 -r 058f2dc7226e scripts/general/integral2.m
--- a/scripts/general/integral2.m	Tue Jul 12 22:19:54 2016 -0700
+++ b/scripts/general/integral2.m	Sun Oct 08 01:42:53 2017 -0400
@@ -1,221 +1,420 @@
-## Copyright (C) 2017 Nicholas Jankowski
-##
-## 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  {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb})
-## @deftypefnx {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{prop}, @var{val}, @dots{})
-##
-## Numerically evaluate the two-dimensional integral of @var{f} using adaptive
-## quadrature over the two-dimensional domain defined by @var{xa}, @var{xb},
-## @var{ya}, @var{yb} (scalars may be finite or infinite).
-##
-## @code{integral2} is a wrapper for @code{dblquad} intended to provide
-## @sc{matlab} compatibility.  More control of the numerical integration may
-## be achievable by calling the various quadrature functions directly.
-##
-## @var{f} is a function handle, inline function, or string containing the name
-## of the function to evaluate.  The function @var{f} must be of the form
-## @math{z = f(x,y)} where @var{x} is a vector and @var{y} is a scalar.  It
-## should return a vector of the same length and orientation as @var{x}.
-##
-## Additional optional parameters can be specified using
-## @qcode{"@var{property}", @var{value}} pairs.  Valid properties are:
-##
-## @table @code
-## @item AbsTol
-## Define the absolute error tolerance for the quadrature.  The default
-## value is 1e-10 (1e-5 for single).
-##
-## @item RelTol
-## Define the relative error tolerance for the quadrature.  The default
-## value is 1e-6 (1e-4 for single).
-## @end table
-##
-## Adaptive quadrature is used to minimize the estimate of error until the
-## following is satisfied:
-## @tex
-## $$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
-## @end tex
-## @ifnottex
-##
-## @example
-## @group
-##   @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|).
-## @end group
-## @end example
-##
-## @end ifnottex
-##
-## Known @sc{matlab} incompatibilities:
-##
-## @enumerate
-## @item
-## If tolerances are left unspecified, and any integration limits or waypoints
-## are of type @code{single}, then Octave's integral functions automatically
-## reduce the default absolute and relative error tolerances as specified
-## above.  If tighter tolerances are desired they must be specified.
-## @sc{matlab} leaves the tighter tolerances appropriate for @code{double}
-## inputs in place regardless of the class of the integration limits.
-#
-## @item
-## @code{integral2} currently only integrates functions over rectangular
-## domains.  Implementing @var{ya} and @var{yb} as functions of @var{x} is a
-## planned future improvement.
-##
-## @item
-## The @qcode{"Method"} property is not yet implemented in Octave due to the
-## lack of a @qcode{"tiled"} integrator implementation.  All integrals are
-## evaluated using an equivalent of the @qcode{"iterated"} method.
-## @end enumerate
-##
-## @seealso{integral, integral3, quad, quadgk, quadv, quadl, quadcc, trapz,
-##          dblquad, triplequad}
-## @end deftypefn
-
-function q = integral2 (f, xa, xb, ya, yb, varargin)
-  ## FIXME: It is possible that a non-rectangular domain could be handled by
-  ##        overlaying the integrand with a boolean mask function such that
-  ##        the integration occurs over a rectangle, but regions outside the
-  ##        desired domain contribute zero to the integral.  This may be an
-  ##        inefficient but acceptable hack to get around the rectangular
-  ##        domain limit without having to rewrite the integrating function.
-
-  ## FIXME: Implement "Method" property to let the user select between iterated
-  ##        and tiled integration.  Tiled integration follows the method of
-  ##        Matlab's quad2d function, currently unimplemented in Octave.
-  ##        Should probably just wait for a quad2d implementation to point the
-  ##        integral2 wrapper to, instead of trying to re-create it here.  The
-  ##        following can be added to the help docstring once it is functional:
-  ## @item Method
-  ## Specifies the two dimensional integration method to be used, with valid
-  ## options being @var{"auto"}, @var{"tiled"}, or @var{"iterated"}.
-  ## @code{integral} will use @var{"auto"} by default, where it will usually
-  ## choose @var{"tiled"} unless any of the integration limits are infinite.
-
-  if (nargin < 5 || (mod (nargin, 2) == 0))
-    print_usage ();
-  endif
-
-  if (! is_function_handle (f))
-    print_usage ();
-  endif
-
-  if (! (isscalar (xa) && isscalar (xb) && isscalar (ya)) && isscalar (yb))
-    print_usage ();
-  endif
-
-  ## Check for single or double limits to set appropriate default tolerance.
-  issingle = isa ([xa, xb, ya, yb], "single");
-
-  ## Set defaults, update with any specified parameters.
-  if (issingle)
-    abstol = 1e-5;
-    reltol = 1e-4;
-  else
-    abstol = 1e-10;
-    reltol = 1e-6;
-  endif
-
-  if (nargin == 5)
-    ## Pass the simplest case directly to integrator.
-    ## Let quadcc function handle input checks.
-    q = dblquad (f, xa, xb, ya, yb, [abstol, reltol], @quadcc);
-
-  else
-    ## Parse options to determine how to call integrator.
-    intmethod = [];
-
-    idx = 1;
-    while (idx < nargin - 5)
-      prop = varargin{idx++};
-      if (! ischar (prop))
-        error ("integral2: property PROP must be a string");
-      endif
-
-      switch (tolower (prop))
-        case "abstol"
-          abstol = varargin{idx++};
-          if (! (isnumeric (abstol) && isscalar (abstol) && abstol >= 0))
-            error ("integral2: AbsTol value must be a numeric scalar >= 0");
-          endif
-
-        case "reltol"
-          reltol = varargin{idx++};
-          if (! (isnumeric (reltol) && isscalar (reltol) && reltol >= 0))
-            error ("integral2: RelTol value must be a numeric scalar >= 0");
-          endif
-
-        case "method"
-          intmethod = varargin{idx++};
-          warning (["integral2: Only 'iterated' method implemented.  ", ...
-                    "Method property ignored."]);
-        otherwise
-          error ("integral2: unknown property '%s'", prop);
-      endswitch
-    endwhile
-
-    q = dblquad (f, xa, xb, ya, yb, [abstol, reltol], @quadcc);
-
-  endif
-
-endfunction
-
-
-%!test
-%! f = @(x, y) x .* y;
-%! assert (integral2 (f, 0, 1, 0, 1), 0.25, 1e-10);
-
-%!test
-%! f = @(x, y) 9 * x.^2 + 15 * y.^2;
-%! assert (integral2 (f, 0, 5, -5, 0, "AbsTol", 1e-9), 5000, 1e-9);
-%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6), 5000, -1e-6);
-%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "AbsTol", 1e-9),
-%!         5000, 1e-9);
-
-## tests from dblquad
-%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-7),
-%!        2*log (2), 1e-7)
-%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6),
-%!        2*log (2), -1e-6)
-%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-8,
-%!        "RelTol", 1e-6), 2*log (2), -1e-6)
-%!assert (integral2 (@(x, y) exp (-x.^2 - y.^2) , -1, 1, -1, 1),
-%!        pi * erf (1).^2, 1e-10)
-
-%!assert (integral2 (@plus, 1, 2, 3, 4), 5, 1e-10)
-
-## Test input validation
-%!error integral2
-%!error integral2 (0, 1 ,2 ,3 ,4)
-%!error integral2 (@plus)
-%!error integral2 (@plus, 1)
-%!error integral2 (@plus, 1, 2)
-%!error integral2 (@plus, 1, 2, 3)
-%!error integral2 (@plus, 1, 2, 3, [4 5])
-%!error integral2 (@plus, 1, 2, 3, "test")
-%!error integral2 (@plus, 1, 2, 3, 4, "foo")
-%!error <unknown property 'foo'>  integral2 (@plus, 1, 2, 3, 4, "foo", "bar")
-%!error <property PROP must be a string> integral2 (@plus, 1, 2, 3, 4, 99, "bar")
-%!error <AbsTol value must be a numeric> integral2 (@plus, 1, 2, 3, 4, "AbsTol", "foo")
-%!error <AbsTol value must be a .* scalar> integral2 (@plus, 1, 2, 3, 4, "AbsTol", [1, 2])
-%!error <AbsTol value must be.* .= 0> integral2 (@plus, 1, 2, 3, 4, "AbsTol", -1)
-%!error <RelTol value must be a numeric> integral2 (@plus, 1, 2, 3, 4, "RelTol", "foo")
-%!error <RelTol value must be a .* scalar> integral2 (@plus, 1, 2, 3, 4, "RelTol", [1, 2])
-%!error <RelTol value must be.* .= 0> integral2 (@plus, 1, 2, 3, 4, "RelTol", -1)
-%!warning <Only 'iterated' method implemented>
-%! q = integral2 (@plus, 0, 1, 0, 1, "Method", "tiled");
+
+## Copyright (C) 2017 Nicholas Jankowski, 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  {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb})
+## @deftypefnx {} {@var{q} =} integral2 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{prop}, @var{val}, @dots{})
+## @deftypefnx {} {[@var{q}, @var{err}] =} integral2 (@dots{})
+##
+## Numerically evaluate the two-dimensional integral of @var{f} using adaptive
+## quadrature over the two-dimensional domain defined by @var{xa}, @var{xb},
+## @var{ya}, @var{yb} (scalars may be finite or infinite). Additionally, 
+## @var{ya} and @var{yb} may be scalar functions of @var{x}, allowing for the
+## integration over non-rectangular domains.
+##
+## @var{f} is a function handle, inline function, or string containing the name
+## of the function to evaluate.  The function @var{f} must be of the form
+## @math{z = f(x,y)} where @var{x} is a vector and @var{y} is a scalar.  It
+## should return a vector of the same length and orientation as @var{x}.
+##
+## Additional optional parameters can be specified using
+## @qcode{"@var{property}", @var{value}} pairs.  Valid properties are:
+##
+## @table @code
+## @item Method
+## Specifies the two dimensional integration method to be used, with valid
+## options being @var{"auto"}, @var{"tiled"}, or @var{"iterated"}.
+## @code{integral} will use @var{"auto"} by default, where it will usually
+## choose @var{"tiled"} unless any of the integration limits are infinite.
+##
+## @item Vectorized
+## Option to enable/disable vectorized integration. False forces octave to use 
+## only scalar inputs when calling the integrand, which enables integrands 
+## @math{f(x,y)} that have not been vectorized and only accept @var{x} and 
+## @var{y} as scalars. Default value is @code{true}.
+##
+## @item AbsTol
+## Define the absolute error tolerance for the quadrature.  The default
+## value is 1e-10 (1e-5 for single).
+##
+## @item RelTol
+## Define the relative error tolerance for the quadrature.  The default
+## value is 1e-6 (1e-4 for single).
+## @end table
+##
+## Adaptive quadrature is used to minimize the estimate of error until the
+## following is satisfied:
+## @tex
+## $$error \leq \max \left( AbsTol, RelTol\cdot\vert q\vert \right)$$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##   @var{error} <= max (@var{AbsTol}, @var{RelTol}*|@var{q}|).
+## @end group
+## @end example
+##
+## @end ifnottex
+##
+## @var{err} is an approximate bound on the error in the integral
+## @code{abs (@var{q} - @var{I})}, where @var{I} is the exact value of the
+## integral.
+##
+## Known @sc{matlab} incompatibilities:
+##
+## @enumerate
+## @item
+## If tolerances are left unspecified, and any integration limits
+## are of type @code{single}, then Octave's integral functions automatically
+## reduce the default absolute and relative error tolerances as specified
+## above.  If tighter tolerances are desired they must be specified.
+## @sc{matlab} leaves the tighter tolerances appropriate for @code{double}
+## inputs in place regardless of the class of the integration limits.
+## @end enumerate
+##
+##
+## Reference: @nospell{L.F. Shampine},
+## @cite{"@sc{matlab} program for quadrature in 2D"}, Applied Mathematics and 
+## Computation, pp. 266--274, Vol 1, 2008.
+##
+## @seealso{integral, integral3, quad, quadgk, quadv, quadl, quadcc, trapz,
+##          dblquad, triplequad}
+## @end deftypefn
+
+function [q, err] = integral2 (f, xa, xb, ya, yb, varargin)
+  
+  if (nargin < 5 || (mod (nargin, 2) == 0))
+    print_usage ();
+  endif
+
+  if (! is_function_handle (f))
+    print_usage ();
+  endif
+
+  if (! (isscalar (xa) && isscalar (xb)))
+    print_usage ();
+  endif
+
+  ## Check for single or double limits to set appropriate default tolerance.
+  issingle = isa ([xa, xb], "single");
+  issingle = issingle || ((! is_function_handle (ya)) && isa (ya, "single"));
+  issingle = issingle || ((! is_function_handle (yb)) && isa (yb, "single"));
+
+  ## Set defaults, update with any specified parameters.
+  if issingle
+    abstol = 1e-5;
+    reltol = 1e-4;
+  else
+    abstol = 1e-10;
+    reltol = 1e-6;
+  endif
+  
+  method = "auto";
+  
+  idx = 1;
+  while (idx < nargin - 5)
+    prop = varargin{idx++};
+    if (! ischar (prop))
+      error ("integral2: property PROP must be a string");
+    endif
+
+    switch (tolower (prop))
+      case "abstol"
+        abstol = varargin{idx++};
+        if (! ((isnumeric (abstol)) && (isscalar (abstol)) && (abstol >= 0)))
+          error ("integral2: AbsTol value must be a numeric scalar >= 0");
+        endif
+
+      case "reltol"
+        reltol = varargin{idx++};
+        if (! ((isnumeric (reltol)) && (isscalar (reltol)) && (reltol >= 0)))
+          error ("integral2: RelTol value must be a numeric scalar >= 0");
+        endif
+
+      case "method"
+        method = tolower (varargin{idx++});
+        if (! any (strcmp (method, {"auto", "iterated", "tiled"})))
+          error ("integral2 : method '%s' unrecognized", method)
+        endif
+      case "vectorized"
+        # option to allow unvectorized functions to be used
+        vectorized = varargin{idx++};
+        if (! islogical (vectorized))
+          error ("integral2: 'vectorized' must be a logical value");
+        elseif (! vectorized)
+		  f = @(x, y) arrayfun(f, x, y);
+		endif      
+      
+      otherwise
+        error ("integral2: unknown property '%s'", prop);
+    endswitch
+  endwhile
+    
+  if strcmp (method, "auto")
+    if ((isinf (xa)) || (isinf (xb)) || ...
+         ((! is_function_handle(ya)) && (isinf (ya))) || ...
+         ((! is_function_handle(yb)) && (isinf (yb))))
+      method = "iterated";
+    else
+      method = "tiled";
+    endif
+  endif
+
+  # check upper and lower bounds of y
+  if (! is_function_handle (ya))
+    if isscalar (ya)
+      ya = @(x) ya * ones (rows (x), columns (x));
+    else
+      error ("integral2: 'ya' must be a constant or a (vectorized) function.")
+    endif
+  endif
+  if (! is_function_handle (yb))
+    if isscalar (yb)
+      yb = @(x) yb * ones (rows (x), columns (x));
+    else
+      error ("integral2: 'ya' must be a constant or a (vectorized) function.")
+    endif
+  endif
+  
+  if strcmp(method, "iterated") 
+    q = outer_iterated (f, xa, xb, ya, yb, abstol, reltol);
+    
+    if nargout == 2
+      warning("integral2: 'iterated' method can not return estimated error");
+      err = 0;
+    endif
+    
+  else
+    [q, err] = tiled_integral2(f, xa, xb, ya, yb, abstol, reltol);
+
+  endif  
+endfunction
+
+function q = outer_iterated (f, xa, xb, ya, yb, abstol, reltol)   
+    inner = @inner_iterated;
+    q = feval (@quadcc, @(x) inner (x, f, ya, yb, abstol, reltol), xa, xb, [abstol, reltol]);
+endfunction
+
+function q = inner_iterated (x, f, ya, yb, abstol, reltol)
+    q = zeros (size (x));
+    for i = 1 : length (x)
+      q(i) = feval (@quadcc, @(y) f(x(i), y), ya(x(i)), yb(x(i)), [abstol, reltol]);
+    endfor
+endfunction
+
+function [q, qerr] = tiled_integral2(f, xa, xb, ya, yb, abstol, reltol)
+    maxiter = 5000;
+    iter = 0;  
+    qaccept = 0;
+    qerraccept = 0;
+    area = (xb - xa);
+    
+    # Initialize tile list
+    tilelist(1) = struct("xa", xa, "xb", xb, "ya", 0, "yb", 1, "q", 0, "qerr", Inf);
+    
+    while (length(tilelist) > 0 && iter++ < maxiter)
+      # Get tile with the largest error
+      [~, indx] = max(tilelist.qerr);
+      tile = tilelist(indx);
+      tilelist(indx) = [];
+      
+      # Subdivide the tile into 4 subtiles
+      tiles(4) = struct("xa", tile.xa, "xb", tile.xa + tile.xb / 2, 
+                        "ya", tile.ya, "yb", tile.ya + tile.yb / 2,
+                        "q", 0, "qerr", 0);
+      tiles(3) = struct("xa", tile.xa, "xb", tile.xa + tile.xb / 2, 
+                        "ya", tile.ya + tile.yb / 2, "yb", tile.yb,
+                        "q", 0, "qerr", 0);
+      tiles(2) = struct("xa", tile.xa + tile.xb / 2, "xb", tile.xb, 
+                        "ya", tile.ya, "yb", tile.ya + tile.yb / 2,
+                        "q", 0, "qerr", 0);
+      tiles(1) = struct("xa", tile.xa + tile.xb / 2, "xb", tile.xb, 
+                        "ya", tile.ya + tile.yb / 2, "yb", tile.yb,
+                        "q", 0, "qerr", 0);
+
+      # Perform the quadrature of 4 subtiles
+      for i = 1:4
+        [tiles(i).q, tiles(i).qerr] = tensorproduct (f, ya, yb, tiles(i));
+      endfor
+ 
+      q = qaccept + sum([[tilelist.q], tiles.q]);
+      qerr = qerraccept + sum([[tilelist.qerr], tiles.qerr]);
+      tol = max (abstol, reltol .* abs (q));
+      localtol = tol * ([tile.xb] - [tile.xa]) * ([tile.yb] - [tile.ya]) / area;
+      
+      # If global tolerance is met, return
+      if (qerr < tol)
+        break;
+      endif
+  
+      # Accept the tiles meeting the tolerance, and add the others back to
+      # the list of tiles to treat
+      indx = find([tiles.qerr] < localtol);
+      qaccept += sum([tiles(indx).q]);
+      qerraccept += sum([tiles(indx).qerr]);
+      tiles(indx) = [];
+      tilelist = [tilelist, tiles];
+    endwhile
+    
+    if (iter >= maxiter)
+      if (qerr > max (abstol, reltol .* abs (q)))
+        warning ("integral2 : Maximum number of sub-tiles without convergence.")
+      else
+        warning("integral2: Maximum number of sub-tiles, maybe low accuracy.")
+      endif
+    endif   
+endfunction
+
+function [q, qerr] = tensorproduct (f, ya, yb, tile)
+  # The Shampine TwoD paper proposes using a G3,K7 rule in a tensor product
+  # I couldn't find a tabulated abscissas et weighrs of a G3,K7 rule publically
+  # available, so I'm using a G7,G15 rule from Octave's implementation of quadgk 
+  
+  persistent abscissa = [-0.9914553711208126e+00, -0.9491079123427585e+00, ...
+                         -0.8648644233597691e+00, -0.7415311855993944e+00, ...
+                         -0.5860872354676911e+00, -0.4058451513773972e+00, ...
+                         -0.2077849550078985e+00,  0.0000000000000000e+00, ...
+                          0.2077849550078985e+00,  0.4058451513773972e+00, ...
+                          0.5860872354676911e+00,  0.7415311855993944e+00, ...
+                          0.8648644233597691e+00,  0.9491079123427585e+00, ...
+                          0.9914553711208126e+00];
+
+  persistent weights15 = ...
+            [0.2293532201052922e-01,  0.6309209262997855e-01, ...
+             0.1047900103222502e+00,  0.1406532597155259e+00, ...
+             0.1690047266392679e+00,  0.1903505780647854e+00, ...
+             0.2044329400752989e+00,  0.2094821410847278e+00, ...
+             0.2044329400752989e+00,  0.1903505780647854e+00, ...
+             0.1690047266392679e+00,  0.1406532597155259e+00, ...
+             0.1047900103222502e+00,  0.6309209262997855e-01, ...
+             0.2293532201052922e-01];
+
+  persistent weights7  = [0.0, ...
+             0.1294849661688697e+00, 0.0, ...
+             0.2797053914892767e+00, 0.0, ...
+             0.3818300505051889e+00, 0.0, ...
+             0.4179591836734694e+00, 0.0, ...
+             0.3818300505051889e+00, 0.0, ...
+             0.2797053914892767e+00, 0.0, ...
+             0.1294849661688697e+00, 0.0];
+             
+  xaa = tile.xa;
+  xbb = tile.xb;
+  yaa = tile.ya;
+  ybb = tile.yb;
+             
+  xcenter = (xbb + xaa) / 2;
+  xhalfwidth =  (xbb - xaa) / 2;  
+  x = xhalfwidth * abscissa + xcenter;
+  xx = ones(15, 1) * x;
+  ylower = ya(x);
+  yupper = yb(x);
+  yhalfwidth = (yupper - ylower) .* (ybb - yaa) ./ 2;
+  ycenter = ylower + (yupper - ylower) .* (ybb + yaa) ./ 2;
+  yy = (abscissa' * yhalfwidth + ones(15, 1) * ycenter);
+  zz = f (xx, yy);
+  
+  q = weights15 * (weights15 * (ones(15,1) * yhalfwidth .* zz))' .* xhalfwidth;
+  qerr = abs(weights7 * (weights7 * (ones(15,1) * yhalfwidth .* zz))' .* xhalfwidth
+        - q);
+endfunction
+
+## method tests
+%!test
+%! f = @(x, y) x .* y;
+%! assert (integral2 (f, 0, 1, 0, 1), 0.25, 1e-10);
+%! assert (integral2 (f, 0, 1, 0, 1, "method", "tiled"), 0.25, 1e-10);
+%! assert (integral2 (f, 0, 1, 0, 1, "method", "iterated"), 0.25, 1e-10);
+%! assert (integral2 (f, 0, 1, 0, 1, "method", "auto"), 0.25, 1e-10);
+
+## vectorized = false test s
+%!test
+%! f = @(x, y) x * y;
+%! assert (integral2 (f, 0, 1, 0, 1, "vectorized", false), 0.25, 1e-10);
+%!error integral2 (f, 0, 1, 0, 1, "vectorized", true)
+
+## tolerance tests
+%!test
+%! f = @(x, y) 9 * x.^2 + 15 * y.^2;
+%! assert (integral2 (f, 0, 5, -5, 0, "AbsTol", 1e-9), 5000, 1e-9);
+%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6), 5000, -1e-6);
+%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "AbsTol", 1e-9),
+%!         5000, 1e-9);
+
+## tests from dblquad
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-7),
+%!        2*log (2), 1e-6);
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6),
+%!        2*log (2), -1e-6);
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-8,
+%!        "RelTol", 1e-6), 2*log (2), -1e-6);
+%!assert (integral2 (@(x, y) exp (-x.^2 - y.^2) , -1, 1, -1, 1),
+%!        pi * erf (1).^2, 1e-10);
+
+%!assert (integral2 (@plus, 1, 2, 3, 4), 5, 1e-10);
+
+%!assert (integral2 (@(x,y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x), 1, -1e-6);
+
+% method="iterated" 
+%!test
+%! f = @(x, y) x .* y;
+%! assert (integral2 (f, 0, 1, 0, 1, "method", "iterated"), 0.25, 1e-10);
+
+%!test
+%! f = @(x, y) 9 * x.^2 + 15 * y.^2;
+%! assert (integral2 (f, 0, 5, -5, 0, "AbsTol", 1e-9, "method", "iterated"), 
+%!         5000, 1e-9);
+%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "method", "iterated"), 
+%!         5000, -1e-6);
+%! assert (integral2 (f, 0, 5, -5, 0, "RelTol", 1e-6, "method", "iterated",
+%!         "AbsTol", 1e-9),5000, 1e-9);
+
+## tests from dblquad
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-7, "method",
+%!         "iterated"), 2*log (2), 1e-7);
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "RelTol", 1e-6, "method",
+%!         "iterated"), 2*log (2), -1e-6);
+%!assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, "AbsTol", 1e-8,
+%!        "RelTol", 1e-6, "method", "iterated"), 2*log (2), -1e-6);
+%!assert (integral2 (@(x, y) exp (-x.^2 - y.^2) , -1, 1, -1, 1, "method",
+%!          "iterated"), pi * erf (1).^2, 1e-10);
+%!assert (integral2 (@plus, 1, 2, 3, 4, "method", "iterated"), 5, 1e-10);
+%!assert (integral2 (@(x,y) 1 ./ (x + y), 0, 1, 0, @(x) 1 - x, "method", 
+%!         "iterated"), 1, -1e-6);
+
+## Test input validation
+%!error integral2
+%!error integral2 (0, 1 ,2 ,3 ,4)
+%!error integral2 (@plus)
+%!error integral2 (@plus, 1)
+%!error integral2 (@plus, 1, 2)
+%!error integral2 (@plus, 1, 2, 3)
+%!error integral2 (@plus, 1, 2, 3, [4 5])
+%!error integral2 (@plus, 1, 2, 3, "test")
+%!error integral2 (@plus, 1, 2, 3, 4, "foo")
+%!error <unknown property 'foo'>  integral2 (@plus, 1, 2, 3, 4, "foo", "bar")
+%!error <property PROP must be a string> integral2 (@plus, 1, 2, 3, 4, 99, "bar")
+%!error <AbsTol value must be a numeric> integral2 (@plus, 1, 2, 3, 4, "AbsTol", "foo")
+%!error <AbsTol value must be a .* scalar> integral2 (@plus, 1, 2, 3, 4, "AbsTol", [1, 2])
+%!error <AbsTol value must be.* .= 0> integral2 (@plus, 1, 2, 3, 4, "AbsTol", -1)
+%!error <RelTol value must be a numeric> integral2 (@plus, 1, 2, 3, 4, "RelTol", "foo")
+%!error <RelTol value must be a .* scalar> integral2 (@plus, 1, 2, 3, 4, "RelTol", [1, 2])
+%!error <RelTol value must be.* .= 0> integral2 (@plus, 1, 2, 3, 4, "RelTol", -1)
+%!error <method 'bad' unrecognized> integral2 (@plus, 1, 2, 3, 4, "method", "bad")
