# HG changeset patch
# User Nicholas R. Jankowski <jankowskin@asme.org>
# Date 1506015517 14400
#      Thu Sep 21 13:38:37 2017 -0400
# Node ID 948356c5704dd788c8a05c84f502b5bdbc95c569
# Parent  1091931bd63caa2c2eff57ffbe2b40a646dbf809
Add missing functions integral2.m and integral3.m (bug #52074)


* integral2.m: new function file, primarily a wrapper to dblquad
* integral3.m: new function file, primarily a wrapper to triplequad
* /scripts/general/module.mk: add integral2.m and integral3.m to the list

* quad.txi: add doc reference to integral2 and 3 after description of dblquad   and triplequad
* NEWS: announce new function integral
2, integral3
* __unimplemented__.m: remove integral2, integral3 from the list

diff -r 1091931bd63c -r 948356c5704d NEWS
--- a/NEWS	Thu Sep 21 12:45:28 2017 -0400
+++ b/NEWS	Thu Sep 21 13:38:37 2017 -0400
@@ -49,6 +49,8 @@
       hgtransform
       humps
       integral
+	  integral2
+	  integral3
       openvar
       repelem
       rticks
diff -r 1091931bd63c -r 948356c5704d doc/interpreter/quad.txi
--- a/doc/interpreter/quad.txi	Thu Sep 21 12:45:28 2017 -0400
+++ b/doc/interpreter/quad.txi	Thu Sep 21 13:38:37 2017 -0400
@@ -318,6 +318,7 @@
 @end group
 @end example
 
+@item dblquad, triplequad
 The above process can be simplified with the @code{dblquad} and
 @code{triplequad} functions for integrals over two and three
 variables.  For example:
@@ -333,6 +334,14 @@
 
 @DOCSTRING(triplequad)
 
+@item integral2, integral3
+@code{integral2} and @code{integral3} are also provided as partial 
+compatibility wrappers to @{dblquad} and @{triplequad}:
+
+@DOCSTRING(integral2)
+
+@DOCSTRING(integral3)
+
 The above mentioned approach works, but is fairly slow, and that problem
 increases exponentially with the dimensionality of the integral.  Another
 possible solution is to use Orthogonal Collocation as described in the
diff -r 1091931bd63c -r 948356c5704d scripts/general/integral2.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/integral2.m	Thu Sep 21 13:38:37 2017 -0400
@@ -0,0 +1,199 @@
+## 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 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
+## absolute tolerance is 1e-10 (1e-5 for single).
+## @end table
+##
+## Known Matlab incompatibilities:
+## @table @code
+## @item
+## 1. @code{integral2} currently only functions over rectangular domains. 
+## Implementing @var{ya} and @var{yb} as functions of @var{x} is a planned 
+## future improvement.
+## @item
+## 2. A @var{'Method'} property is not yet implemented in Octave due to the lack
+## of a 'tiled' integration implementation. All integrals are evaluated using an
+## equivalent of the 'iterated' method. 
+## @item
+## 3. The underlying 2d integrator only accepts an Absolute Tolerance. As such 
+## it is not possible to specify @var{'RelTol'}. A default Relative Tolerance
+## value of 1e-6 (1e-4 for single) is used unless Relative Tolerance testing is
+## disabled by specifying @var{'RelTol'} as @var{'off'}.
+## 
+## @end table
+##
+## @seealso{quad, quadgk, quadv, quadl, quadcc, trapz, intergral, 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 than 
+  ##        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 recreate 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. 
+
+  ## FIXME: implement 'reltol' property once there is a good way to pass this 
+  ##        value to the underlying integrator. The following can be added to
+  ##        the help docstring once it is functional:
+  ## @item RelTol
+  ## Define the relative error tolerance for the quadrature.  The default
+  ## relative tolerance is 1e-6 (1e-4 for single).
+
+  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 so can set appropriate default tolerance
+  issingle = (isa (xa, "single") || isa (xb, "single") || ... 
+              isa (ya, "single") || 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
+
+  intmethod = [];
+  integfunc = @quadgk;
+
+  ## check optional parameters to adjust defaults
+  
+  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 "reltol"
+        reltol = varargin{idx++};
+        if (! ischar (reltol))
+          warning("integral2: RelTol cannot currently be changed from ", ...
+                  " the default. It can be disabled by specifying ", ...\
+                  "'RelTol', 'off'.");
+        endif
+      case "abstol"
+        abstol = varargin{idx++};
+        if (! isscalar(abstol))
+          error("integral2: AbsTol value must be a scalar")
+        endif
+        
+      case "method"
+        intmethod = varargin{idx++};
+        warning (["integral2: alternate integration methods not yet ", ...
+                  "implemented. Method property ignored."]);
+      otherwise
+        error ("integral2: unknown property '%s'", prop);
+    endswitch
+  endwhile
+    
+  if (ischar (reltol))
+    if (strcmp (tolower (reltol), "off"))
+      #if don't want reltol defined, set integrator to quadcc
+      integfunc = @quadcc;
+    else
+      error ("integral2: unknown RelTol value '%s'", reltol);
+    endif
+  endif
+
+  q = dblquad (f, xa, xb, ya, yb, abstol, integfunc);
+
+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', 'off'), 5000, 1e-10);
+
+## tests from dblquad
+%! assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, 'AbsTol', 1e-7), 
+%!         2*log (2), 1e-9);
+%! assert (integral2 (@(x, y) 1 ./ (x+y), 0, 1, 0, 1, 'RelTol', 'off'), 
+%!         2*log (2), 1e-10);
+%! 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 (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 scalar> integral2 (@plus, 1, 2, 3, 4, 'AbsTol', 'foo')
+%!error <unknown RelTol value 'foo'> integral2 (@plus, 1, 2, 3, 4, 'RelTol', 'foo')
\ No newline at end of file
diff -r 1091931bd63c -r 948356c5704d scripts/general/integral3.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/integral3.m	Thu Sep 21 13:38:37 2017 -0400
@@ -0,0 +1,202 @@
+## 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} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb})
+## @deftypefnx {} {@var{q} =} integral3 (@var{f}, @var{xa}, @var{xb}, @var{ya}, @var{yb}, @var{za}, @var{zb}, @var{prop}, @var{val}, @dots{})
+##
+## Numerically evaluate the three dimensional integral of @var{f} using adaptive 
+## quadrature over the three-dimensional domain defined by @var{xa}, @var{xb}, 
+## @var{ya}, @var{yb}, @var{za}, @var{zb} (scalars may be finite or infinite).
+##
+## @code{integral3} is a wrapper for @code{triplequad} intended to provide 
+## 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{w = f(x,y,z)} where either @var{x} or @var{y} is a vector and the 
+## remaining inputs are scalars. @var{f} should return a vector of the same 
+## length and orientation as the vector input @var{x} or @var{y}.
+##
+## 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
+## absolute tolerance is 1e-10 (1e-5 for single).
+## @end table
+##
+## Known Matlab incompatibilities:
+## @table @code
+## @item
+## 1. @code{integral3} currently only functions over a rectangular volume. 
+## Implementing @var{ya} and @var{yb} as functions of @var{x}, and @var{za} and
+## @var{zb} as functions of (@var{x,y})is a planned future improvement.
+##
+## @item
+## 2. A @var{'Method'} property is not yet implemented in Octave due to the lack
+## of a 'tiled' integration implementation. All integrals are evaluated using an
+## equivalent of the 'iterated' method. 
+##
+## @item
+## 3. The underlying 3d integrator only accepts an Absolute Tolerance. As such 
+## it is not possible to specify @var{'RelTol'}. A default Relative Tolerance
+## value of 1e-6 (1e-4 for single) is used unless Relative Tolerance testing is
+## disabled by specifying @var{'RelTol'} as @var{'off'}.
+## 
+## @end table
+##
+## @seealso{quad, quadgk, quadv, quadl, quadcc, trapz, intergral, dblquad,
+##  triplequad}
+## @end deftypefn
+
+function q = integral3 (f, xa, xb, ya, yb, za, zb, varargin)
+  ## FIXME: it is possible that a non-rectangular domain could be handled by 
+  ##        overlaying the integrand with a boolean mask function such than 
+  ##        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 
+  ##        integral3 wrapper to instead of trying to recreate 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. 
+
+  ## FIXME: implement 'reltol' property once there is a good way to pass this 
+  ##        value to the underlying integrator. The following can be added to
+  ##        the help docstring once it is functional:
+  ## @item RelTol
+  ## Define the relative error tolerance for the quadrature.  The default
+  ## relative tolerance is 1e-6 (1e-4 for single).
+
+  if (nargin < 7 || (mod (nargin, 2) == 0))
+    print_usage ();
+  endif
+  
+  if (! is_function_handle (f))
+    print_usage ();
+  endif
+   
+  if ((! isscalar (xa)) || (! isscalar (xb)) ...
+      || (! isscalar (ya)) || (! isscalar (yb)) ...
+      || (! isscalar (za)) || (! isscalar (zb)))
+    print_usage ();
+  endif
+  
+  #check for single or double limits so can set appropriate default tolerance
+  issingle = (isa (xa, "single") || isa (xb, "single") || ... 
+              isa (ya, "single") || isa (yb, "single") || ...
+              isa (za, "single") || isa (zb, "single"));
+  
+  ## Set defaults, update with any specified parameters.
+  if issingle
+    abstol = 1e-5;
+    reltol = 1e-4;  
+  else
+    abstol = 1e-10;
+    reltol = 1e-6;
+  endif
+
+  intmethod = [];
+  integfunc = @quadgk;
+
+  ## check optional parameters to adjust defaults
+  
+  idx = 1;
+  while (idx < nargin - 7)
+    prop = varargin{idx++};
+    if (! ischar (prop))
+      error ("integral3: property PROP must be a string");
+    endif
+
+    switch (tolower (prop))
+      case "reltol"
+        reltol = varargin{idx++};
+        if (! ischar (reltol))
+          warning("integral3: RelTol cannot currently be changed from ", ...
+                  " the default. It can be disabled by specifying ", ...\
+                  "'RelTol', 'off'.");
+        endif
+      case "abstol"
+        abstol = varargin{idx++};
+        if (! isscalar(abstol))
+          error("integral3: AbsTol value must be a scalar")
+        endif
+        
+      case "method"
+        intmethod = varargin{idx++};
+        warning (["integral3: alternate integration methods not yet ", ...
+                  "implemented. Method property ignored."]);
+      otherwise
+        error ("integral3: unknown property '%s'", prop);
+    endswitch
+  endwhile
+    
+  if (ischar (reltol))
+    if (strcmp (tolower (reltol), "off"))
+      #if don't want reltol defined, set integrator to quadcc
+      integfunc = @quadcc;
+    else
+      error ("integral3: unknown RelTol value '%s'", reltol);
+    endif
+  endif
+
+  q = triplequad (f, xa, xb, ya, yb, za, zb, abstol, integfunc);
+
+endfunction
+
+%!test
+%! f = @(x, y, z) x.*y.*z;
+%! assert (integral3 (f, 0, 1, 0, 1, 0, 1), 0.125, 1e-10);
+
+%!test
+%! f = @(x,y,z) y.*sin(x)+z.*cos(x);
+%! assert (integral3 (f, 0, pi, 0, 1, -1, 1), 2, 1e-10);
+%! assert (integral3 (f, 0, pi, 0, 1, -1, 1, 'RelTol', 'off'), 2, 1e-10);
+
+## tests from triplequad
+%! assert (integral3 (@(x,y,z) exp (-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1,
+%!          1), pi^(3/2) * erf (1).^3, 1e-10);
+%! assert (integral3 (@(x,y,z) exp (-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1,
+%!           1, "RelTol", "off"), pi^(3/2) * erf (1).^3, 1e-10);
+
+## Test input validation
+%!error integral3 (0, 1 ,2 ,3 ,4, 5, 6)
+%!error integral3 (@plus)
+%!error integral3 (@plus, 1)
+%!error integral3 (@plus, 1, 2)
+%!error integral3 (@plus, 1, 2, 3)
+%!error integral3 (@plus, 1, 2, 3, 4)
+%!error integral3 (@plus, 1, 2, 3, 4, 5)
+%!error integral3 (@plus, 1, 2, 3, 4, 5, [6 7])
+%!error integral3 (@plus, 1, 2, 3, 4, 5, 'test')
+%!error integral3 (@plus, 1, 2, 3, 4, 5, 6, 'foo')
+%!error <unknown property 'foo'>  integral3 (@plus, 1, 2, 3, 4, 5, 6, 'foo', 'bar')
+%!error <property PROP must be a string> integral3 (@plus, 1, 2, 3, 4, 5, 6, 99, 'bar')
+%!error <AbsTol value must be a scalar> integral3 (@plus, 1, 2, 3, 4, 5, 6, 'AbsTol', 'foo')
+%!error <unknown RelTol value 'foo'> integral3 (@plus, 1, 2, 3, 4, 5, 6, 'RelTol', 'foo')
\ No newline at end of file
diff -r 1091931bd63c -r 948356c5704d scripts/general/module.mk
--- a/scripts/general/module.mk	Thu Sep 21 12:45:28 2017 -0400
+++ b/scripts/general/module.mk	Thu Sep 21 13:38:37 2017 -0400
@@ -42,6 +42,8 @@
   %reldir%/inputParser.m \
   %reldir%/int2str.m \
   %reldir%/integral.m \
+  %reldir%/integral2.m \
+  %reldir%/integral3.m \
   %reldir%/interp1.m \
   %reldir%/interp2.m \
   %reldir%/interp3.m \
diff -r 1091931bd63c -r 948356c5704d scripts/help/__unimplemented__.m
--- a/scripts/help/__unimplemented__.m	Thu Sep 21 12:45:28 2017 -0400
+++ b/scripts/help/__unimplemented__.m	Thu Sep 21 13:38:37 2017 -0400
@@ -59,12 +59,6 @@
       txt = ["griddedInterpolant is not implemented.  ", ...
              "Consider using griddata."];
 
-    case "integral2"
-      txt = ["integral2 is not implemented.  Consider using dblquad."];
-
-    case "integral3"
-      txt = ["integral3 is not implemented.  Consider using triplequad"];
-
     case "linprog"
       txt = ["Octave does not currently provide linprog.  ", ...
              "Linear programming problems may be solved using @code{glpk}.  ", ...
