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

ENH: Root finding #417

Merged
merged 3 commits into from Jul 25, 2018

Conversation

Projects
None yet
6 participants
@chrishyland

chrishyland commented Jul 4, 2018

Hi @jstac, in this PR @spvdchachan and I added in the Jitted Newton methods for root finding, In particular, the newton method and the newton_secant method. Additionally, we wrote some test cases too.

As several people mentioned in #416 , extra information about convergence, iterations, etc. is crucial.
Therefore, we provided this extra information but instead used the namedtuple collection for fields to be accessed easier through names rather than numerical indexes in normal tuples.

We followed the source code directly from Scipy but there were several technical issues with Numba.

  • Optional function args: Numba only supports optional parameters of primitive types only. Also, it is not guaranteed to be 100% stable as well (explicitly mentioned in the doc). This means it is not possible at the moment to support all three variants of newton's method within one jitted function. i.e. we cannot make fprime and fprime2 optional. Currently, we've only implemented 2 of the functions.
  • Warnings: Cannot issue a warning message through warnings. If this is crucial, we have to find a workaround.

@jstac

spvdchachan added some commits Jul 4, 2018

Add jitted newton methods for root finding
Based on Scipy's newton methods, we add two variants of newton root
finding - Newton Raphson and secant method. All methods are jitted
through Numba with nopython mode. Relevant information (apart from the
root) is also returned in the result as python namedtuple since Numba
supports it.
@coveralls

This comment has been minimized.

Show comment
Hide comment
@coveralls

coveralls Jul 4, 2018

Coverage Status

Coverage increased (+0.04%) to 95.142% when pulling 876c797 on chrishyland:root-finding into d6130c1 on QuantEcon:master.

coveralls commented Jul 4, 2018

Coverage Status

Coverage increased (+0.04%) to 95.142% when pulling 876c797 on chrishyland:root-finding into d6130c1 on QuantEcon:master.

@jstac

This comment has been minimized.

Show comment
Hide comment
@jstac

jstac Jul 5, 2018

Contributor

Great, thanks. Nice looking code.

@mmcky using namedtuple leads to an extra dependency. Are you OK with that? I am in favor of using it.

If we do use named tuples then let's also change our existing maximization function accordingly.

I can't think of a better approach to accommodating the fact that only only primitive arguments can be optional. @oyamad Do you have any suggestions on how to better structure this code?

This is a trivial point, but I thought I read somewhere that having multiple return statements within one function is frowned on these days in software circles. It that not so? In any case you might also save a few bytes by replacing your multiple return statements with break.

Contributor

jstac commented Jul 5, 2018

Great, thanks. Nice looking code.

@mmcky using namedtuple leads to an extra dependency. Are you OK with that? I am in favor of using it.

If we do use named tuples then let's also change our existing maximization function accordingly.

I can't think of a better approach to accommodating the fact that only only primitive arguments can be optional. @oyamad Do you have any suggestions on how to better structure this code?

This is a trivial point, but I thought I read somewhere that having multiple return statements within one function is frowned on these days in software circles. It that not so? In any case you might also save a few bytes by replacing your multiple return statements with break.

@oyamad

This comment has been minimized.

Show comment
Hide comment
@oyamad

oyamad Jul 5, 2018

Member

@jstac I am not sure if I really understand the issue, but the following generates different code for different signatures (which behaves like optional arguments from the user's perspective):

from numba import jit, generated_jit, types

@generated_jit(nopython=True)
def test(x, f=None, f2=None):
    if isinstance(f, types.Callable):
        if isinstance(f2, types.Callable):
            def return_f2(x, f, f2):
                return f2(x)
            return return_f2
        # else
        def return_f(x, f, f2):
            return f(x)
        return return_f
    
    #else
    def return_0(x, f, f2):
        return 0
    return return_0
@jit(nopython=True)
def f(x):
    return x

@jit(nopython=True)
def f2(x):
    return 2*x
test(10)
0
test(10, f)
10
test(10, f, f2)
20

But, in my opinion, this would unnecessarily complicate the structure of the code and increase the maintenance costs. Having different names for different sets of arguments (newton for (func, x0, fprime), newton_secant (func, x0), and newton_??? for (func, x0, fprime, fprime2)) appears fine to me.

  • I don't think issuing a warning message is crucial. Just returning the information about whether the routine is converged is fine.
Member

oyamad commented Jul 5, 2018

@jstac I am not sure if I really understand the issue, but the following generates different code for different signatures (which behaves like optional arguments from the user's perspective):

from numba import jit, generated_jit, types

@generated_jit(nopython=True)
def test(x, f=None, f2=None):
    if isinstance(f, types.Callable):
        if isinstance(f2, types.Callable):
            def return_f2(x, f, f2):
                return f2(x)
            return return_f2
        # else
        def return_f(x, f, f2):
            return f(x)
        return return_f
    
    #else
    def return_0(x, f, f2):
        return 0
    return return_0
@jit(nopython=True)
def f(x):
    return x

@jit(nopython=True)
def f2(x):
    return 2*x
test(10)
0
test(10, f)
10
test(10, f, f2)
20

But, in my opinion, this would unnecessarily complicate the structure of the code and increase the maintenance costs. Having different names for different sets of arguments (newton for (func, x0, fprime), newton_secant (func, x0), and newton_??? for (func, x0, fprime, fprime2)) appears fine to me.

  • I don't think issuing a warning message is crucial. Just returning the information about whether the routine is converged is fine.
Add newton halley method and modify returning
The third and last newton method for root finding is now added with its
test cases. Each newton method is modified to include only a single
return statement.
@spvdchachan

This comment has been minimized.

Show comment
Hide comment
@spvdchachan

spvdchachan Jul 8, 2018

Contributor

Thank you for the input.
@oyamad I've added another newton method called newton_halley which can take in fprime2.

@jstac I've removed multiple returns in each method as requested.

Contributor

spvdchachan commented Jul 8, 2018

Thank you for the input.
@oyamad I've added another newton method called newton_halley which can take in fprime2.

@jstac I've removed multiple returns in each method as requested.

@jstac

This comment has been minimized.

Show comment
Hide comment
@jstac

jstac Jul 23, 2018

Contributor

@oyamad @mmcky This looks good to me. I would be happy to see it merged.

I would also like to see a robust root finding method such as brent or bisection. But perhaps this can be done separately.

Contributor

jstac commented Jul 23, 2018

@oyamad @mmcky This looks good to me. I would be happy to see it merged.

I would also like to see a robust root finding method such as brent or bisection. But perhaps this can be done separately.

@mmcky

This comment has been minimized.

Show comment
Hide comment
@mmcky

mmcky Jul 24, 2018

Contributor

thanks for your review @jstac and @oyamad. Thanks @chrishyland and @spvdchachan for this PR.
@jstac I have no problem with adding additional dependencies such as named tuples

Contributor

mmcky commented Jul 24, 2018

thanks for your review @jstac and @oyamad. Thanks @chrishyland and @spvdchachan for this PR.
@jstac I have no problem with adding additional dependencies such as named tuples

@mmcky mmcky added the ready label Jul 24, 2018

@mmcky mmcky changed the title from Root finding to ENH: Root finding Jul 24, 2018

@mmcky mmcky added the enhancement label Jul 24, 2018

@mmcky

This comment has been minimized.

Show comment
Hide comment
@mmcky

mmcky Jul 25, 2018

Contributor

thanks again for everyone's input. I will merge this PR now.

Contributor

mmcky commented Jul 25, 2018

thanks again for everyone's input. I will merge this PR now.

@mmcky mmcky merged commit 1155267 into QuantEcon:master Jul 25, 2018

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment