## 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).
##
## @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).
##
## @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 disable vectorized integration, forcing octave to use only scalar 
## inputs when calling the integrand.
## @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

## The "titled" method of this code is a conversion of the integral2.R function
## from the R package pracma and inherents its GPLv3 license

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
  
  ## FIXME
  ## Shampine included the possibility of non vectorized functions. The 
  ## code below supports this option, but it don't seem to be used in the
  ## matlab integral2 function. Decide how to treat it; either the option
  ## or document it
  vectorized = true;
  
  singular = false;
  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 = varargin{idx++};
               
      case "vectorized"
        # Undocumented option to allow unvectorized functions to be treated
        vectorized = varargin{idx++};
        if !islogical (vectorized)
          error ("integral2: 'vectorized' must be a logical value");
        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
  elseif !(strcmp(method, "iterated") || strcmp(method, "tiled"))
    error ("integral2 : method '%s' unrecognized", method)
  endif

  if ! vectorized
    f = @(x, y) arrayfun(f, x, y);
  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

## 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")

% method="auto" which defaults to 'tiled'
%!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-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" which defaults to 'tiled'
%!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-6)
%!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)