# HG changeset patch
# User Nicholas R. Jankowski <jankowskin@asme.org>
# Date 1504625121 14400
#      Tue Sep 05 11:25:21 2017 -0400
# Node ID c49bf57fb68117e44a8e8b20ba5b2a9cb12f6fda
# Parent  c2fbb639f970c9b032bfae9f89fe08afc664e692
Add missing function integral.m (bug #42037)

* integral.m: new function file, primarily a wrapper to quadgk or quadv
* /scripts/general/module.mk: add integral.m to the list
* quad.txi: add doc reference to integral after description of quad functions
* NEWS: annouce new function integral
* __unimplemented__.m: remove integral from the list

diff -r c2fbb639f970 -r c49bf57fb681 NEWS
--- a/NEWS	Sun Sep 17 21:45:01 2017 +0200
+++ b/NEWS	Tue Sep 05 11:25:21 2017 -0400
@@ -48,6 +48,7 @@
       gsvd
       hgtransform
       humps
+	  integral
       openvar
       repelem
 
diff -r c2fbb639f970 -r c49bf57fb681 doc/interpreter/quad.txi
--- a/doc/interpreter/quad.txi	Sun Sep 17 21:45:01 2017 +0200
+++ b/doc/interpreter/quad.txi	Tue Sep 05 11:25:21 2017 -0400
@@ -32,14 +32,14 @@
 @node Functions of One Variable
 @section Functions of One Variable
 
-Octave supports five different algorithms for computing the integral
+Octave supports five different adaptive quadrature algorithms for computing
+the integral
 @tex
 $$
  \int_a^b f(x) d x
 $$
 @end tex
-of a function @math{f} over the interval from @math{a} to @math{b}.
-These are
+of a function @math{f} over the interval from @math{a} to @math{b}. These are
 
 @table @code
 @item quad
@@ -57,6 +57,12 @@
 @item quadcc
 Numerical integration using adaptive @nospell{Clenshaw-Curtis} rules.
 
+In addition, the following functions are also provided:
+
+@item integral
+A compatibility wrapper function that will choose between @code(quadv) and
+@code(quadgk) depending on the integrand and options chosen.
+
 @item trapz, cumtrapz
 Numerical integration of data using the trapezoidal method.
 @end table
@@ -165,6 +171,8 @@
 
 @DOCSTRING(quadcc)
 
+@DOCSTRING(integral)
+
 Sometimes one does not have the function, but only the raw (x, y) points from
 which to perform an integration.  This can occur when collecting data in an
 experiment.  The @code{trapz} function can integrate these values as shown in
diff -r c2fbb639f970 -r c49bf57fb681 scripts/general/integral.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/integral.m	Tue Sep 05 11:25:21 2017 -0400
@@ -0,0 +1,209 @@
+## Copyright (C) 2017 Nicholas Jankowski
+##
+## This program 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.
+##
+## This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {} {@var{q} =} integral (@var{f}, @var{a}, @var{b})
+## @deftypefnx {} {@var{q} =} integral (@var{f}, @var{a}, @var{b}, @var{prop1}, @var{val1}, @var{prop2}, @var{val2}, @dots{})
+##
+## Numerically evaluate the integral of @var{f} from @var{a} to @var{b}
+## using adaptive quadrature.
+##
+## @code{integral} is a wrapper for @code{quadgk} (for scalar integrands) and
+## @code{quadv} (for array-vauled integrands) 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 vectorized and
+## return a vector of output values when given a vector of input values.
+##
+## @var{a} and @var{b} are the lower and upper limits of integration.  Either
+## or both limits may be infinite or contain weak end singularities. If either
+## or both limits are complex, @code{integral} will perform a straight line path
+## integral.  Alternatively, a complex domain path can be specified using the
+## "Waypoints" option (see below).
+##
+## Additional optional parameters can be specified using {"@var{prop}",
+## @var{val}} 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).
+##
+## @item RelTol
+## Define the relative error tolerance for the quadrature.  The default
+## relative tolerance is 1e-6 (1e-4 for single).
+##
+## @item Waypoints
+## Specifies points to be used in defining subintervals of the quadrature
+## algorithm, or if @var{a}, @var{b}, or @var{waypoints} are complex then
+## the quadrature is calculated as a contour integral along a piecewise
+## continuous path. For more detail see @code{quadgk}.
+##
+## @item ArrayValued
+## @code{integral} expects @var{f} to returns a scalar value unless
+## @var(arrayvalued) is specified as @code{true}. This will cause
+## @code{integral} to perform the integration over the entire array and
+## return @var{q} with the same dimensions as returned by @var{f}.
+## @end table
+##
+## *Note: As a consequence of using @code{quadgk} and @code{quadv}, certain
+## option combinations are currently unsupported. "ArrayValued" cannot be
+## combined with "RelTol" or "Waypoints". This is a known incompatibility with
+## Matlab.
+##
+## @seealso{quad, quadgk, quadv, quadl, quadcc, trapz, dblquad, triplequad}
+## @end deftypefn
+
+function Q = integral (F, A, B, varargin)
+
+  if (nargin < 3 || (mod (nargin, 2) == 0))
+    print_usage ();
+  endif
+
+  ## Pass simplest case directly to genereral integrator. minimize overhead
+  ## and let quad function handle input checks on function and limits.
+
+  if (nargin == 3)
+    Q = quadgk (F, A, B);
+
+  else
+    ## Parse options to determine how to call integrator
+    abstol = [];
+    reltol = [];
+    waypoints = [];
+    arrayvalued = false;
+
+    idx = 1;
+    while (idx < nargin - 3)
+      prop = varargin{idx++};
+      if (! ischar (prop))
+        error ("integral: property PROP must be a string");
+      endif
+      
+      switch (tolower (prop))
+        case "reltol"
+          reltol = varargin{idx++};
+        case "abstol"
+          abstol = varargin{idx++};
+        case "waypoints"
+          waypoints = varargin{idx++}(:);
+        case "arrayvalued"
+          arrayvalued = varargin{idx++};
+        otherwise
+          error ("integral: unknown property '%s'", prop);
+      endswitch
+    endwhile
+
+    if (arrayvalued)
+      ## FIXME: replace warning with arrayfun(?) call to quadgk
+      if (! isempty (waypoints))
+        warning(["integral: array-valued quadrature routine currently ", ...
+                 "unable to handle WayPoints. WayPoints are ignored."]);
+      endif
+
+      ## FIXME: remove warning once we have reltol compatible arrayval'd
+      ## quad or arrayfun call to quadgk.
+      if (! isempty (reltol))
+        warning(["integral: array-valued quadrature only accepts AbsTol. ", ...
+                 "RelTol ignored."]);
+      endif
+
+      ## quadv accepts empty abstol input.
+      Q = quadv (F, A, B, abstol);
+
+    else
+      ## quadgk will accept empty values of optional parameters
+      Q = quadgk (F, A, B, "AbsTol", abstol, "RelTol", reltol, ...
+                  "WayPoints", waypoints);
+
+    endif
+  endif
+endfunction
+
+%!demo
+%! a = 0:2;
+%! b = 2*pi:pi:4*pi;
+%! integral (@sin, a, b)
+
+## Matlab compatibility tests
+%!test
+%! f = @(x) exp(-x.^2).*log(x).^2;
+%! emgamma = 0.57721566490153286060651209008240243104215933593992;
+%! ansexact = (sqrt(pi)*(8*log(2)^2+8*emgamma*log(2)+pi^2+2*emgamma^2))/16;
+%! assert (integral (f, 0, inf), ansexact, 1e-6);
+%! assert (integral (f, 0, inf, "RelTol", 1e-12), ansexact, 1e-12);
+
+%!test ##with parameter:
+%! f = @(x,c) 1./(x.^3-2*x-c);
+%! assert (integral (@(x) f(x,5), 0, 2), -0.4605015338467329, 1e-12);
+
+%!test ##with tolerances
+%! f = @(x) log(x);
+%! assert (integral (@(x) f(x), 0, 1, "AbsTol", 1e-6), -1, 1e-6);
+
+%!test ##waypoints
+%! f = @(x) 1./(2.*x-1);
+%! assert (integral (f, 0, 0,"Waypoints",[1+1i, 1-1i]), -pi*1i, 1e-12);
+
+%!test ##test array function
+%! f = @(x) sin ((1:5)*x);
+%! assert (integral(f, 0, 1, "ArrayValued", true),1./[1:5]-cos(1:5)./[1:5],1e-6);
+
+%!test
+%! f = @(x) x .^ 5 .* exp (-x) .* sin (x);
+%! assert (integral (f, 0, inf, "RelTol", 1e-8, "AbsTol", 1e-12), -15, 2e-12);
+
+##tests from quadgk
+%!test
+%!assert (integral (@sin,-pi,pi), 0, 1e-6)
+%!assert (integral (inline ("sin"),-pi,pi), 0, 1e-6)
+%!assert (integral ("sin",-pi,pi), 0, 1e-6)
+%!assert (integral (@sin,-pi,pi), 0, 1e-6)
+
+%!assert (integral (@sin,-pi,0), -2, 1e-6)
+%!assert (integral (@sin,0,pi), 2, 1e-6)
+%!assert (integral (@(x) 1./sqrt (x),0,1), 2, 1e-6)
+%!assert (integral (@(x) abs (1 - x.^2),0,2, "Waypoints", 1), 2, 1e-6)
+%!assert (integral (@(x) 1./(sqrt (x) .* (x+1)),0,Inf), pi, 1e-6)
+%!assert (integral (@(z) log (z),1+1i,1+1i, "WayPoints", [1-1i, -1,-1i, -1+1i]), -pi * 1i, 1e-6)
+%!assert (integral (@(x) exp (-x .^ 2),-Inf,Inf), sqrt (pi), 1e-6)
+%!assert (integral (@(x) exp (-x .^ 2),-Inf,0), sqrt (pi)/2, 1e-6)
+
+
+## tests from quadv
+## Test vector-valued functions
+%!assert (integral (@(x) [(sin (x)), (sin (2 * x))], 0, pi, "ArrayValued", true)
+%!         , [2, 0], 1e-5)
+
+## Test matrix-valued functions
+%!test
+%! warning ("off", "Octave:divide-by-zero", "local");
+%! assert (integral (@(x) [x, x, x; x, 1./sqrt(x), x; x, x, x], 0, 1,
+%!         "ArrayValued", true), [0.5, 0.5, 0.5; 0.5, 2, 0.5; 0.5, 0.5, 0.5],
+%!         1e-5);
+
+## Test input validation
+%!error integral (@sin)
+%!error integral (@sin, 0)
+%!error integral (@sin, 0, 1, 1e-6, true, 4)
+%!error integral (@sin, 0, 1, "DummyArg")
+%!error <property PROP must be a string> integral (@sin, 0, 1, "AbsTol", 1e-6, 2, 3)
+%!error <unknown property 'foo'> integral (@sin, 0, 1, "foo", 3)
+%!error integral (@sin, 0, 1, "AbsTol", ones (2,2))
+%!error integral (@sin, 0, 1, "AbsTol", -1)
+%!error integral (@sin, 0, 1, "RelTol", ones (2,2))
+%!error integral (@sin, 0, 1, "RelTol", -1)
\ No newline at end of file
diff -r c2fbb639f970 -r c49bf57fb681 scripts/general/module.mk
--- a/scripts/general/module.mk	Sun Sep 17 21:45:01 2017 +0200
+++ b/scripts/general/module.mk	Tue Sep 05 11:25:21 2017 -0400
@@ -41,6 +41,7 @@
   %reldir%/idivide.m \
   %reldir%/inputParser.m \
   %reldir%/int2str.m \
+  %reldir%/integral.m \
   %reldir%/interp1.m \
   %reldir%/interp2.m \
   %reldir%/interp3.m \
diff -r c2fbb639f970 -r c49bf57fb681 scripts/help/__unimplemented__.m
--- a/scripts/help/__unimplemented__.m	Sun Sep 17 21:45:01 2017 +0200
+++ b/scripts/help/__unimplemented__.m	Tue Sep 05 11:25:21 2017 -0400
@@ -59,10 +59,6 @@
       txt = ["griddedInterpolant is not implemented.  ", ...
              "Consider using griddata."];
 
-    case "integral"
-      txt = ["Octave provides many routines for 1-D numerical integration.  ", ...
-             "Consider quadcc, quad, quadv, quadl, quadgk."];
-
     case "integral2"
       txt = ["integral2 is not implemented.  Consider using dblquad."];
 
