New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/issue 210 numerical integration #566

Merged
merged 40 commits into from Feb 7, 2018

Conversation

Projects
None yet
9 participants
@randommm
Member

randommm commented Jun 22, 2017

Submission Checklist

  • Run unit tests: ./runTests.py test/unit/math/prim/arr/functor/integrate_function_test.cpp test/unit/math/rev/arr/functor/integrate_function_test.cpp
  • Run cpplint: make cpplint
  • Declare copyright holder and open-source license: see below

Summary:

Fix issue #210

Intended Effect:

An easier way to do numerical integration on Stan with the potential benefit of allowing the user to define gradient function which might make things faster on some cases.

How to Verify:

There is some unit tests included.

Side Effects:

There's a potential problem: I'm unsure about the relative tolerance: the original implementation of John D. Cook only has absolute tolerance:

            if (errorEstimate < 0.1*targetAbsoluteError)
                break;

What I did was adding a simple relative tolerance measure:

            previousIntegral = integral;
            integral = 0.5*integral + newContribution;
            ...
            if (errorEstimate < 0.1*targetAbsoluteError ||
                fabs(1 - integral / previousIntegral) <
                targetRelativeError)
                break;

But unsure how appropriate and safe this is. For example, should it wait for at least 10 iterations of the algorithm:

            previousIntegral = integral;
            integral = 0.5*integral + newContribution;
            ...
            if (errorEstimate < 0.1*targetAbsoluteError)
                break;

            if (fabs(1 - integral / previousIntegral) <
                targetRelativeError && levels >= 10)
                break;

Documentation:

Included.

It might also be the case to change the name of this function to some thing like tanh_sinh to allow make way for possibly other integration methods in the future? I'm open to suggestion here.

Reviewer Suggestions:

Possibly @bgoodri for the problem described above and maybe @bob-carpenter or @syclik on the C++ code side?

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company):

By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:

randommm added some commits Dec 3, 2015

It seems that boost::bind placeholders are declared in namespace unna…
…med and this is messing up with Stan.

According to the author of boost::bind here: http://boost.2283326.n4.nabble.com/reference-to-1-is-ambiguous-td2589379.html boost::bind is able to use boost::lambda placeholders, so this solves the problem.
Merge branch 'develop' into feature/issue-210-numerical-integration
 Conflicts:
	stan/math.hpp
	stan/math/prim/mat/fun/value_of.hpp

@randommm randommm requested review from syclik, bob-carpenter and bgoodri Jun 22, 2017

@Eh2406

This comment has been minimized.

Eh2406 commented Jun 22, 2017

Note that tanh_sinh does not work well if the integrand violates the assumptions. So big +1 on clearly documenting that we are using tanh_sinh.

@seantalts

This comment has been minimized.

Member

seantalts commented Sep 21, 2017

Jenkins, retest this please.

@randommm

This comment has been minimized.

Member

randommm commented Sep 25, 2017

OK, @bob-carpenter it's ready for review.

Note: I'm not sure how to deal with indentation and formatting of the lines with de_integrator(std::bind<double>, I tried to make it in way that is easier to read.

@bob-carpenter

This only needs a few minor things:

  • a bit of doc
  • change an include format
  • rearrange a loop
    I'm very excited about getting this into the language since a lot of people have been asking for it.
*/
template <typename F, typename G, typename T_param>
inline typename scalar_type<T_param>::type
integrate_1d_tsc_tscg(const F& f, const G& g, const double a,

This comment has been minimized.

@bob-carpenter

bob-carpenter Sep 26, 2017

Contributor

There is no explanation of what tsc or tscg are. The integration algorithm and conditions on it also need to be described. What kinds of functions can be integrated? Do the bounds need to be finite? Are the absolute and relative error such that both must hold or is one enough?

This comment has been minimized.

@randommm

randommm Oct 24, 2017

Member

I think what cause the confusion was that this was meant to be integrate_1d_tscg, but a search and replace mistake made it integrate_1d_tsc_tscg. Fixed.

* double (double, T_param) where the first argument is one being
* integrated and the second one is either an extra scalar or vector
* being passed to f.
* @param g a functor with signature

This comment has been minimized.

@bob-carpenter

bob-carpenter Sep 26, 2017

Contributor

Why are funtor and gradient being passed separately? I'd think they'd be much easier to compute together.

This comment has been minimized.

@randommm

randommm Oct 1, 2017

Member

You mean requesting F to return a vector with its value and gradients all together? Makes a lot of sense and a lot of problems might have common structures between them which could reused to greatly speed up calculations. I'm a little busy this week, so I'll try to finish this PR next weekend.

This comment has been minimized.

@randommm

randommm Oct 24, 2017

Member

Unfortunatelly I checked that this is not possible, at least not in the current implementation: the function and every single one of its gradients must be integrated separately.

#ifndef STAN_MATH_REV_MAT_FUNCTOR_DEINTEGRATOR_HPP
#define STAN_MATH_REV_MAT_FUNCTOR_DEINTEGRATOR_HPP
#include "de_integrator_constants.hpp"

This comment has been minimized.

@bob-carpenter

bob-carpenter Sep 26, 2017

Contributor

This should follow the pattern of our other includes with <stan/math/rev/mat/functor/de_integrator_constants.hpp>.

namespace math {
/**
* Double Exponential Integrator.

This comment has been minimized.

@bob-carpenter

bob-carpenter Sep 26, 2017

Contributor

Is there a citation for the algorithm being used?

Are there conditions on the function f?

Can the upper and lower limits be infiite?

Must absolute and relative error requirements both be met or is one enough?

return ops_partials.build(val_);
// easy case, here we are calculating a normalizing constant,
// not a normalizing factor, so g doesn't matter at all
} else {

This comment has been minimized.

@bob-carpenter

bob-carpenter Sep 26, 2017

Contributor

You want to reverse the conditions, return, then not indent. So rather than

if (!special_case) {
  A
} else {
  B
}

use

if (special_case)
  B;
A;

It removes the level of nesting and the level of negation. Both make the overall structure of the main algorithm clearer.

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Sep 26, 2017

Just to be clear here, I'm not asking for a long fancy algorithm description. Just a spelled-out name and a link to a citation is enough.

@syclik

This comment has been minimized.

Member

syclik commented Oct 22, 2017

@bob-carpenter, @randommm. Status on this pull request?

@randommm

This comment has been minimized.

Member

randommm commented Oct 22, 2017

@syclik

This comment has been minimized.

Member

syclik commented Oct 22, 2017

No need to apologize! I'm looking forward to pull request being in the code base.

@seantalts

This comment has been minimized.

Member

seantalts commented Nov 1, 2017

(Please ignore the 3rd failing check - still working out the kinks)

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Jan 30, 2018

@yizhang-cae Should we close this pending your addition of an integrator?

@randommm

This comment has been minimized.

Member

randommm commented Jan 30, 2018

@yizhang-cae

This comment has been minimized.

Contributor

yizhang-cae commented Jan 30, 2018

@bob-carpenter I don't CVODES integrator addresses the issue here. With CVODES ODE solvers cannot do improper integrals such as in (-infinity, +infinity).

@yizhang-cae

This comment has been minimized.

Contributor

yizhang-cae commented Jan 30, 2018

@randommm I've added another ODE solver(Adams-Moulton scheme from CVODES) for the need of high-order integrator. But I don't think it addresses the issue here. IIUC here we need 3 things that are different from current ODE solver(but in the background the implementation may be based on ODE solvers):

  • Different functor signature, with the functor simply being a univariate integrand
  • Able to do improper integral
  • Integral limits as parameters
@betanalpha

This comment has been minimized.

Contributor

betanalpha commented Jan 30, 2018

@randommm

This comment has been minimized.

Member

randommm commented Jan 30, 2018

@betanalpha I addressed all the requested changes by @bob-carpenter, so I thought that all that was left was a final review and merge, and then working the parser out on stan-dev/stan repo.

I believe item 3 of @yizhang-cae said could be released in a future version of this integrator (it would be backward compatible), I don't know if item 2 is "compatible" with the algorithm, and unsure about item 1.

@betanalpha

This comment has been minimized.

Contributor

betanalpha commented Jan 30, 2018

@yizhang-cae

This comment has been minimized.

Contributor

yizhang-cae commented Jan 30, 2018

yes, this is how quadpack does the infty integel.

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Jan 30, 2018

@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Feb 7, 2018

Ack, I see this is waiting for my final review. If I miss stuff like this, please let me know!

@bob-carpenter bob-carpenter merged commit a4c9203 into develop Feb 7, 2018

2 checks passed

continuous-integration/jenkins/pr-merge This commit looks good
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
@bob-carpenter

This comment has been minimized.

Contributor

bob-carpenter commented Feb 7, 2018

OK, now I need to add the Stan issue to match before this is properly exposed to the language.

@seantalts seantalts deleted the feature/issue-210-numerical-integration branch Feb 7, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment