bugGNU Octave - Bugs: bug #42037, Missing function integral

 
 

bug #42037: Missing function integral

Submitter:  Nir Krakauer <nir_krakauer>
Submitted:  Thu 03 Apr 2014 01:01:43 PM UTC
   
 
Category:  Octave Function Severity:  1 - Wish
Priority:  5 - Normal Item Group:  Feature Request
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 17 Oct 2017 08:32:54 PM UTC, comment #29: 

Indeed, I had forgotten about this closed report.  I reviewed and checked in your changes here (http://hg.savannah.gnu.org/hgweb/octave/rev/47dd094a6239).

Some of the BIST tests were now failing so I re-did them and grabbed the BIST tests from quadcc, those with "WayPoints" from quadgk, and kept the Array-Valued tests from quadv.

I changed to using a simple if statement rather than using ifelse since it is more efficient unless the computed branches are particularly simple, like a constant.


Rik <rik5>
Group administrator
Tue 17 Oct 2017 06:28:47 PM UTC, comment #28: 

oops, always a second too late. I see that patch version had the wrong bug number. fixed here in v3

(file #42185)

Nicholas Jankowski <nrjank>
Group Member
Tue 17 Oct 2017 06:16:59 PM UTC, comment #27: 

tweaked version of the last patch that preserves type single output.

(file #42184)

Nicholas Jankowski <nrjank>
Group Member
Wed 27 Sep 2017 09:01:00 PM UTC, comment #26: 

patch updating integral to make more use of new quadcc attached. again, I'm not set up to compile and actually use the new quadcc until it appears in a released windows build, so you may want to check the tests before pushing this.

(file #41903)

Nicholas Jankowski <nrjank>
Group Member
Tue 26 Sep 2017 04:34:39 PM UTC, comment #25: 

ok, now that the integral2 and 3 patch is pushed, I'll tweak integral to use quadcc sometime soon.

Nicholas Jankowski <nrjank>
Group Member
Mon 25 Sep 2017 08:13:14 PM UTC, comment #24: 

I think so.  Not only is quadcc faster than quadgk, because it is a compiled function rather than an m-file, but it can readily handle non-smooth integrands with weird discontinuities including NaNs.  If a user knows something about the specifics of the integrand, such as it being smooth and continuous, then they can choose the integrator that best suits the problem.  But if you don't know anything about it, and only have a handle to the function, you need an integrator that can handle strange behavior.

Rik <rik5>
Group administrator
Mon 25 Sep 2017 07:35:28 PM UTC, comment #23: 

now that quadcc takes both abs and reltol, would it be preferred to mod this function to use quadcc over quadgk if the requested parameters permit? (i.e., use quadcc unless waypoints (then use quadgk) or arrayvalued(then use quadv)...

Nicholas Jankowski <nrjank>
Group Member
Thu 21 Sep 2017 02:46:52 PM UTC, comment #22: 

remembering now. opened a new bug report for integral2 and 3. will detail any issues over there. See Bug #52074

Nicholas Jankowski <nrjank>
Group Member
Wed 20 Sep 2017 05:43:07 PM UTC, comment #21: 

re: dblquad - I initially thought so too when I made this one, but ran into a few problems and stopped. I don't recall what.  I'll review and maybe open a new bug report once i have an idea.

Nicholas Jankowski <nrjank>
Group Member
Wed 20 Sep 2017 05:30:20 PM UTC, comment #20: 

I checked in a slightly different change here (http://hg.savannah.gnu.org/hgweb/octave/rev/6570fdb7d3a0).  A lot of users will be scientists, not mathematicians, so rather than give them vague advice as Matlab does, "Single precision integrations may require larger tolerances", I chose to use the sensible defaults that quadgk has.  They can always try again if their integration fails, but this will allow most integrands to complete successfully.

Also, I redid the precisions in the BIST tests.  The tolerance argument of the assert() function is an absolute tolerance when TOL is positive and a relative tolerance when TOL is negative.  I used this fact to match a BIST test where AbsTol is specified to a positive TOL for assert() and vice versa if RelTol was specified.

I think enough progress has been made to close this report.

Incidentally, integral2 can probably be implemented similarly as a wrapper around dblquad.

Rik <rik5>
Group administrator
Wed 20 Sep 2017 04:27:10 PM UTC, comment #19: 

quadv has a default abstol of 1e-6 which matches matlab. so for compatibility reasons I think quadv should stay as is.

Simply passing abstol = 1e-10 to that integral test has it pass, and best as I can tell versus matlab 2017a this is compatible behavior.

Also, the different error values for double and single are unique to quadgk best as I can tell, and matches matlab quadgk default tolerances according to the help.

There is no specification of different tolerance limits based on class in the help for integral. I had just kept that text there from quadgk, but integral documentation simply states tha abstol default is 1e-10, and RelTol is 1e-6, and "Single precision integrations may require larger tolerances."

My thoughts on the simplest thing to do that also maintains compatibility is just to pass 1e-10 and 1e-6 down to the quad functions (throwing a reltol unimplemented warning for quadv until that gets fixed), and remove the 'single' distinction in the help text.

The attached update_integral.diff patch (builds on the checked in version) makes that change to default reltol and abstol and the help text. Also, added tests for checking single vs double input and output to maintain compatibility.  last, I switched the matrix-valued test from 1/sqrt(x) to exp(x) to get rid of the singularity warnings during test. a couple minor spacing fixes too.

Also changed my Octave editor setting from CRLF (default on windows) to LF.


(file #41845)

Nicholas Jankowski <nrjank>
Group Member
Wed 20 Sep 2017 12:46:10 AM UTC, comment #18: 

I made some other changes including using the default Copyright notice for Octave, using spaces between the function name and the opening parenthesis, using lower case for variable names such as q, f, a, b, and using only one '#' for single-line comments appearing after code.  I checked everything in here (http://hg.savannah.gnu.org/hgweb/octave/rev/2eae2ad53eb9).

I changed one of the array-valued tests to use the default AbsTol of 1e-10 which now causes it to fail.  Although quadgk uses a default of 1e-10, quadv uses a default of 1e-6.

The integral() function could be changed to pass abstol of 1e-10 down to quadv, but I think it would be better to change quadv to simply use these new default values.  For one thing, quadv does not distinguish between double and single values and always uses 1e-6 by default.  That seems wrong.  There is code in the input validation for quadgk which could be copied to quadv to make this work. 


  issingle = (isa (a, "single") || isa (b, "single")
              || isa (waypoints, "single"));

  if (isempty (abstol))
    abstol = ifelse (issingle, 1e-5, 1e-10);
  elseif (! isscalar (abstol) || abstol < 0)
    error ("quadgk: ABSTOL must be a scalar >=0");
  endif



Rik <rik5>
Group administrator
Tue 19 Sep 2017 09:21:11 PM UTC, comment #17: 

Sorry, please find attached a new revised version in which I removed a (failing) demo, leftover from a brave (but aborted) intention to add demos :-).

I also tried to make tests fit into 80 columns and added ";" at the end of each "assert" call.

(file #41838)

Pantxo Diribarne <pantxo>
Group Member
Tue 19 Sep 2017 09:54:36 AM UTC, comment #16: 

Hi,

I reviewed the coding style (see attached) :

  • Not operator: replaced  ~ by !
  • Comments: replaced % by # and removed a few redundant comments
  • Encoding: changed to utf-8 with unix LF
  • Removed the "author" entry at the end of the help text


I also removed "integral" from the list of unimplemented functions.

I did not change the comments about eventually calling quadgk using arrayfun but I really don't see how this would be possible: quadgk expects a function that accepts vector inputs while "arrayvalued" is supposed to be used with functions that accept scalar inputs.

I will push the patch in the coming days unless someone wants to have a last look.

(file #41835)

Pantxo Diribarne <pantxo>
Group Member
Thu 07 Sep 2017 06:03:00 PM UTC, comment #15: 

I have reviewed mainly its functional behavior (checking the BISTs and some corner cases), and less
coding style and performance issues.

Avinoam Kalma <avinoam>
Group Member
Wed 06 Sep 2017 10:37:18 PM UTC, comment #14: 

@Avinoam: I think I saw that you had already given one review.  Is that correct?  That means the second review before committing can be shorter.

Rik <rik5>
Group administrator
Wed 06 Sep 2017 05:13:31 PM UTC, comment #13: 

Thanks,

@Maintainers: Can you review and push this patch?

Avinoam Kalma <avinoam>
Group Member
Tue 05 Sep 2017 03:46:55 PM UTC, comment #12: 

oops.  here's a patch that includes the function.

(file #41734)

Nicholas Jankowski <nrjank>
Group Member
Tue 05 Sep 2017 03:43:26 AM UTC, comment #11: 

Integral.m seems to be missing

Avinoam Kalma <avinoam>
Group Member
Tue 05 Sep 2017 03:41:01 AM UTC, comment #10: 

v2 patch attached that fixes the style issues.

(file #41729)

Nicholas Jankowski <nrjank>
Group Member
Thu 31 Aug 2017 10:52:21 PM UTC, comment #9: 

Will check and reupload with style cleanups sometime this weekend

The MATLAB and math works comments just note a section with example problems to solve for compatibility, taken from public facing help and documentation. One of them is integral of log(x)... No copying of code or similar. Could remove or come up with alternate expressions but don't see any issue testing it with the same expressions.

I would prefer any fixes to the inf error be in quadgk and keep this a wrapper. If numerical fixes appear in here it could be a mess to keep track if later if changes happen in quadgk.

If we can make a test do the same with quadgk, and IR the other quads, it might be worth a separate bug report.

Nicholas Jankowski <nrjank>
Group Member
Thu 31 Aug 2017 03:32:19 PM UTC, comment #8: 

Thanks for implementing integral.

I have checked integral function on Win-7

All test pass:


>> test integral
PASSES 31 out of 31 tests


Coding style:
1. In the test code, some of the keywords are given with double quotes ("RelTol") and some with single quotes ('AbsTol')
2. Missing some spaces after commas.

Using Matlab examples -  Is it legal? What about copyright issues?

Problem with Inf:


>> fun = @(x) 1./x
fun =

@(x) 1 ./ x

>> a = integral (fun, 0, 1)
warning: quadgk: non finite integrand encountered
warning: called from
    quadgk at line 341 column 7
    integral at line 87 column 7
warning: quadgk: Error tolerance not met.  Estimated error 2.41923
warning: called from
    quadgk at line 397 column 5
    integral at line 87 column 7
a =  38.5083785045586


Matlab result:


>> fun = @(x) 1./x

fun =

  function_handle with value:

    @(x)1./x

>> a = integral (fun, 0, 1)
Warning: Infinite or Not-a-Number value encountered.
> In integralCalc/iterateScalarValued (line 349)
  In integralCalc/vadapt (line 132)
  In integralCalc (line 75)
  In integral (line 88)

a =

   Inf


38 is not Inf. Probably the problem is in quadgk and not in integral

What about a simple demo, like:


fun = @(x) cos (x);
integ = @(x) integral (fun, 0, x);
z = [0:0.01:2*pi];
w = arrayfun (integ, z);
plot (z, cos(z), z, w);
title ("integral demo");
h = legend ("f(x) = cos (x)", "integral of f(x)");
legend (h, "location", "southwest");
axis ("tight");
grid on



Avinoam Kalma <avinoam>
Group Member
Tue 21 Mar 2017 08:07:00 PM UTC, comment #7: 

Here's a first attempt at a wrapper function for integral. Function and patch attached.

For compatibility integral needs to handle only 4 options: RelTol, AbsTol, Waypoints, and Arrayvalued. 

If option "ArrayValued" is true, it calls quadv.

Otherwise it will call quadgk. one issue with this approach is that you cannot call ArrayValued along with Waypoints or RelTol. (quadv only takes abstol).  at the moment I have it continue with a warning. these could be changed to errors if you think that's a better handling method.

Also, it's possible the combination could be handled by calling arrayfun and quadgk for those cases, but I wasn't able to come up with a form for arrayfun that actually worked.

In any case, this passes all tests and should provide some bit of improved compatibility. haven't looked at whether adding similar wrappers for integral2 and integral3 will be as straightforward.

(file #40052, file #40053)

Nicholas Jankowski <nrjank>
Group Member
Thu 15 Dec 2016 02:34:26 PM UTC, comment #6: 

This issue is still valid in Octave 4.2.0.

(There is still no Octave implementation of integral.)

Hartmut <hardy>
Wed 16 Jul 2014 03:12:51 AM UTC, comment #5: 

Sorry, you're quite right.  I thought quadcc had been extended to include WayPoints.

Rik <rik5>
Group administrator
Wed 16 Jul 2014 02:03:17 AM UTC, comment #4: 


>The "WayPoint" feature is also implemented in quadcc.


I cannot find such description in help and manual concerning quadcc.

So I tried: (octave-3.8.1 on windows 7)


octave> quadgk (@(x) abs (1 - x.^2), 0, 2, "Waypoints", 1)
ans =  2.0000


I simply substitute quadgk by quadcc. The result is


octave> quadcc (@(x) abs (1 - x.^2), 0, 2, "Waypoints", 1)
error: quadcc: tolerance (TOL) must be a single real scalar > 0


It seems that "WayPoints" feature of quadcc seems not to be implemented in the same way for quadgk.

I pull the code of quadcc.cc in the source tree and search by the keyword 'Waypoints' but nothing is found.

I also the error messages in quadcc.cc,

quadcc: integrand F must return a single, real-valued vector"
quadcc: integrand F must return a single, real-valued vector of the same size as the input

MATLAB function 'integral' supports complex integral so that it can be used only real case.

My observation, only quadgk can handle complex integral in octave.

Tatsuro MATSUOKA <tmacchant>
Tue 15 Jul 2014 03:02:35 PM UTC, comment #3: 

The "WayPoint" feature is also implemented in quadcc which is a good default integration algorithm because it handles discontinuities, Infinities, WayPoints, as well as simple continuous functions.

Rik <rik5>
Group administrator
Tue 15 Jul 2014 07:09:32 AM UTC, comment #2: 

I have looked into octave code for quadl, quadv, and quadgk

The 'Waypoints' feature is implemented in quadgk.

Sorry for my poor knowledge for this matter.

Tatsuro MATSUOKA <tmacchant>
Tue 15 Jul 2014 04:52:25 AM UTC, comment #1: 

The similar discussion is also posted in octave-maintainers ML.
http://octave.1599824.n4.nabble.com/quad-will-be-removed-in-a-future-release-of-MATLAB-td4665467.html

MATLAB integral seem to be not just renaming the function from quad to integral.


q = integral(fun,xmin,xmax,Name,Value)


is new feature of MATLAB integral.  In particular, when Name is 'Waypoints', the implementation seems not to be straightforward.

The 'Waypoints' supports complex integral via complex integral path.
http://www.mathworks.co.jp/jp/help/matlab/ref/integral.html?lang=en#btdd9y9

Tatsuro MATSUOKA <tmacchant>
Thu 03 Apr 2014 01:01:43 PM UTC, original submission:  

Matlab is deprecating the quad function for numerical integration in favor of integral: http://www.mathworks.com/help/matlab/ref/integral.html

I suppose we should add it for compatibility.

Nir Krakauer <nir_krakauer>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #42185:  integral_use_new_quadcc_V3.diff added by nrjank (6KiB - application/octet-stream)
file #42184:  integral_use_new_quadcc_V2.diff added by nrjank (6KiB - application/octet-stream - uses quadcc and preserves single type)
file #41903:  integral_use_new_quadcc.diff added by nrjank (5KiB - application/octet-stream)
file #41845:  update_integral.diff added by nrjank (5KiB - application/octet-stream)
file #41838:  addintegral5.patch added by pantxo (11KiB - text/x-patch)
file #41835:  addintegral4.patch added by pantxo (11KiB - text/x-patch)
file #41734:  add_integral_v3.diff added by nrjank (11KiB - application/octet-stream)
file #41729:  add_integral_v2.diff added by nrjank (2KiB - application/octet-stream)
file #40052:  add_integral.diff added by nrjank (10KiB - application/octet-stream)
file #40053:  integral.m added by nrjank (8KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by pantxo (Updated the item)
  • -email is unavailable- added by avinoam (Updated the item)
  • -email is unavailable- added by nrjank (Updated the item)
  • -email is unavailable- added by hardy (Posted a comment)
  • -email is unavailable- added by mtmiller (Updated the item)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by tmacchant (Posted a comment)
  • -email is unavailable- added by nir_krakauer (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 19 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-10-17 nrjank Attached File- Added integral_use_new_quadcc_V3.diff, #42185
    2017-10-17 nrjank Attached File- Added integral_use_new_quadcc_V2.diff, #42184
    2017-09-27 nrjank Attached File- Added integral_use_new_quadcc.diff, #41903
    2017-09-20 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2017-09-20 nrjank Attached File- Added update_integral.diff, #41845
    2017-09-20 rik5 StatusPatch Reviewed In Progress
    2017-09-19 pantxo Attached File- Added addintegral5.patch, #41838
    2017-09-19 pantxo StatusPatch Submitted Patch Reviewed
    2017-09-19 pantxo Attached File- Added addintegral4.patch, #41835
    2017-09-05 nrjank Attached File- Added add_integral_v3.diff, #41734
    2017-09-05 nrjank Attached File- Added add_integral_v2.diff, #41729
    2017-03-24 avinoam StatusNone Patch Submitted
    2017-03-21 nrjank Attached File- Added add_integral.diff, #40052
        Attached File- Added integral.m, #40053
    2014-11-15 mtmiller CategoryNone Octave Function
        Severity3 - Normal 1 - Wish
        Release3.8.0 dev
        Operating SystemGNU/Linux Any

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code