
## 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.
## @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 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
## For the "iterated" method, @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.
## @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 "singular"
        # Undocumented option to force weakening of edge singularities, but only
        # supported by the 'tiled' method
        singular = varargin{idx++};
        if !islogical (singular)
          error ("integral2: 'singular' must be a logical value");
        endif
        
      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
  
  if strcmp(method, "iterated") 
    if (! (isscalar (ya) && isscalar (yb)))
      error ("integral2: Non scalar limited currently unsupported with the 'iterated' method")
    endif
  
    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
    if !singular
      [q, err, singular] = shampine2(f, xa, xb, ya, yb, abstol, reltol, 
                  issingle, false);
    endif
    if singular
      [q, err] = shampine2(f, xa, xb, ya, yb, abstol, reltol, 
                  issingle, true);
    endif 
  endif  
endfunction

function q = outer_iterated (f, xa, xb, ya, yb, abstol, reltol)
    # 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
    
    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, err, singular] = shampine2(f, xa, xb, ya, yb, abstol, reltol, 
                      issingle, singular)
    # FIXME
    # Should maxlist and maxiter be configurable ?    
    maxlist = 5000;
    maxiter = 5000;
      
    # check input parameters
    nList = floor (maxlist / 10);

    # check upper and lower bounds of y
    if (is_function_handle(ya))
        phiBvar = ya;
    elseif isscalar (ya)
        phiBvar = @(x) ya * ones(rows(x), columns(x));
    else
        error ("integral2: 'ya' must be a constant or a (vectorized) function.")
    endif

    if (is_function_handle(yb))
        phiTvar = yb;
    elseif isscalar (yb)
        phiTvar = @(x) yb * ones(rows(x), columns(x));
    else
        error ("integral2: 'yb' must be a constant or a (vectorized) function.")
    endif
    
    # check borders and redefine
    if (singular)
        thetaL = 0; 
        thetaR = pi;
        phiB = 0;
        phiT = pi;
    else
        thetaL = xa;
        thetaR = xb;
        phiB = 0; 
        phiT = 1;
    endif
    area = (thetaR - thetaL) * (phiT - phiB);
    
    # initial quadrature
    [Qsub, esub] = tensor(xa, xb, thetaL, thetaR, phiB, phiT, f, phiBvar, 
                      phiTvar, singular);
    q = sum (Qsub);

    # some more parameters and main list
    if issingle
      myeps = eps ("single");
    else
      myeps = eps;
    endif
    tol = 100 * eps * abs(q);
    err_ok = 0;
    adjust = 1;
    mainList = zeros(nList, 7);
    nList = 0;
    iter = 1;

    # save info on rectangles in main list
    [mainList, nList, errbnd, err_ok] = save2list(mainList, nList, Qsub, esub,
                        thetaL, thetaR, phiB, phiT, tol, area, adjust, err_ok, 
                        myeps);

    if (nList == 0 || errbnd <= tol)
      err = errbnd;
      return
    endif

    while true
        [mainList, nList, temp] = nextEntry(mainList, nList);

        qa = temp(1); 
        e = temp(2);
        thetaL = temp(3); 
        thetaR = temp(4);
        phiB = temp(5); 
        phiT = temp(6);

        # Approximate integral over four subrectangles
        [Qsub, esub] = tensor(xa, xb, thetaL, thetaR, phiB, phiT, f, phiBvar,
                        phiTvar, singular);

        newq = sum(Qsub);
        adjust = min(1, abs(qa - newq)/e);
        q = q + (newq - qa);
        tol = max(abstol, reltol * abs(q)) / 8;  
        tol = max(tol, 100 * myeps * abs(q));

        [mainList, nList, errbnd2, err_ok] = save2list(mainList, nList, Qsub, esub,
                            thetaL, thetaR, phiB, phiT, tol, area, adjust, err_ok, 
                            myeps);
                            
        # FIXME
        # Try to detect singularities when we aren't weakening them
        # and return early. We'll recall shampine2 with singularity
        # weakening. It seems to work, but is a bit of an ugly test
        if ! singular
          singular = errbnd < errbnd2;
          if singular 
            err = errbnd2;
            return
          endif
        endif
        errbnd = errbnd2;
                            
        if (nList == 0 || errbnd <= tol)
            break;
        elseif (nList > maxlist || iter++ > maxiter)
            if (errbnd > max(abstol, max(100*myeps, reltol) * abs(q)))
              warning ("integral2 : Maximum number of subintervals without convergence.")
            else
              warning("integral2: Maximum number of subintervals, maybe low accuracy.")
            endif
            break
        endif
    endwhile
    err = errbnd;
endfunction

function [qsub, esub] = tensor (a, b, thetaL, thetaR, phiB, phiT, FUN, phiBvar, 
                        phiTvar, singular = false)
    # Gauss-Kronrod (3,7) pair with degrees of precision 5 and 11
    nodes = [-0.9604912687080202, -0.7745966692414834, -0.4342437493468026, ...
                0, 0.4342437493468026, 0.7745966692414834, 0.9604912687080202]; 
    nnodes = length (nodes);
    onevec = ones (2 * nnodes, 1);
    narray = 0.25 * [nodes, nodes];
    wt3 = [0, 5/9, 0, 8/9, 0, 5/9, 0];
    wt7 = [0.1046562260264672, 0.2684880898683334, 0.4013974147759622, ...
           0.4509165386584744, 0.4013974147759622, 0.2684880898683334, ...
           0.1046562260264672];

    qsub = zeros (4, 1);
    esub = zeros (4, 1);
    
    dtheta = thetaR - thetaL;
    etheta = thetaL + dtheta * [0.25, 0.75];
    theta = dtheta*narray + repmat(etheta, nnodes, 1)(:)';
    
    if singular
        x = 0.5 * (b + a) + 0.5 * (b - a) * cos(theta);
    else
        x = theta;
    endif
    X = onevec * x;

    dphi = phiT - phiB;
    ephi = phiB + dphi * [0.25, 0.75];
    phi  = (dphi * narray + repmat(ephi, nnodes, 1)(:)')';

    top = phiTvar(x);
    bottom = phiBvar(x);
    dydt = top - bottom;
    if singular
        t = 0.5 + 0.5*cos(phi);
    else
        t = phi;
    endif
    Y = onevec * bottom + t * dydt;

    Z = FUN(X, Y);
    
    if singular
        temp = 0.25 .* (b - a) .* sin(phi) * (dydt .* sin(theta));
    else
        temp = onevec * dydt;
    endif
    Z = Z .* temp;
    
    # Tensor product: Gauss 3 and 7 points formulae
    esub(1) = wt3 * (wt3 * Z(1:nnodes,1:nnodes))';
    esub(2) = wt3 * (wt3 * Z(1:nnodes,(nnodes+1):(2*nnodes)))';  
    esub(3) = wt3 * (wt3 * Z((nnodes+1):(2*nnodes),1:nnodes))';
    esub(4) = wt3 * (wt3 * Z((nnodes+1):(2*nnodes),(nnodes+1):(2*nnodes)))';
    esub = (esub ./ 4) * (dtheta ./ 2) * (dphi ./ 2);

    qsub(1) = wt7 * (wt7 * Z(1:nnodes,1:nnodes))';
    qsub(2) = wt7 * (wt7 * Z(1:nnodes,(nnodes+1):(2*nnodes)))';
    qsub(3) = wt7 * (wt7 * Z((nnodes+1):(2*nnodes),1:nnodes))';
    qsub(4) = wt7 * (wt7 * Z((nnodes+1):(2*nnodes),(nnodes+1):(2*nnodes)))';
    qsub = (qsub ./ 4) * (dtheta ./ 2) * (dphi ./ 2);  
    
    esub = abs(esub - qsub);
endfunction

function [mainList, nList, errbnd, err_ok] = save2list (mainList, nList,
                      Qsub, esub, thetaL, thetaR, phiB, phiT,
                      tol, area, adjust, err_ok, myeps)

    dtheta = thetaR - thetaL;
    dphi = phiT - phiB;
    localtol = tol * (dtheta/2) * (dphi/2) / area;
    localtol = max(localtol, 100*myeps*abs(sum(Qsub)));

    adjerr = adjust * esub;
    
    if (nList+4) > size(mainList,1)
      mainList = [mainList; zeros(100, 7)];
    endif
    if (adjerr(1) > localtol)
        nList = nList + 1;
        mainList(nList, :) = [Qsub(1), esub(1), thetaL, thetaL + dtheta/2, ...
                                             phiB ,phiB + dphi/2, adjerr(1)];
    else
      err_ok = err_ok + adjerr(1);
    endif

    if (adjerr(2) > localtol)
      nList = nList + 1;
      mainList(nList, :) = [Qsub(2), esub(2), thetaL + dtheta/2, thetaR, ...
                                           phiB, phiB + dphi/2, adjerr(2)];
    else
      err_ok = err_ok + adjerr(2);
    endif

    if (adjerr(3) > localtol)
      nList = nList + 1;
      mainList(nList, :) = [Qsub(3), esub(3), thetaL, thetaL + dtheta/2, ...
                                           phiB + dphi/2, phiT, adjerr(3)];
    else
      err_ok = err_ok + adjerr(3);
    endif

    if adjerr(4) > localtol
      nList = nList + 1;
      mainList(nList, :) = [Qsub(4), esub(4), thetaL + dtheta/2, thetaR, ...
                                           phiB + dphi/2, phiT, adjerr(4)];
    else
      err_ok = err_ok + adjerr(4);
    endif
    
    errbnd = err_ok + sum(mainList(:, 7));
endfunction

function [mainList, nList, entry] = nextEntry (mainList, nList)
    [~, indx] =  max(abs(mainList(1:nList, 7)));
    entry = mainList (indx, :);
    mainList (indx, :) = [];
    nList = nList - 1;
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-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)

%!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-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)