Jim Martin | 18 Apr 19:43 2015

Improvements (scipy.integrate._ode) -- accessing the dense output of DOPRI5 and DOP853.

** Summary **

The wrappers around the differential equation solvers DOPRI5 and DOP853 in "scipy/integrate/_ode.py" do not currently allow access to their "dense output" option (iout=2).  This option provides the coefficients of an interpolating polynomial for the solution between the last and current integration step.  Dense output interpolation is commonly used for event location and improved plotting.  The "_ode.py" code can be modified to optionally allow access to the dense output of DOPRI5/DOP853.  The proposed changes are here:
with example usage here: 
Existing code using "_ode.py" would not be impacted by the changes.

** Background **

The efficient numerical solution of differential equations requires adaptive stepsizes.  Modern high order Runge-Kutta code --- such as DOPRI5 and DOP853 --- can take large variable stepsizes, so that plots made simply by connecting lines between steps are of poor visual quality.  As discussed by [Shampine, 1997, see Section 5 and Figure 1, and Shampine et al., 2003, pages 54, 55] this issue is addressed in the Matlab solvers by interpolation between integration steps, with a default of 4 output solutions provided per step (but can be modified using the "Refine" option). Interpolation of the solution output between steps is also desirable for event location (see discussion in Shampine et al., 2003, Section 2.3.1).

The Fortran codes DOPRI5/DOP853 have a dense output option (iout=2), as discussed in Hairer et al.'s book [Hairer, 1993].  This dense output option is illustrated in the example usage code for these routines:  dr_dopri5.f and dr_dop853.f (see http://www.unige.ch/~hairer/software.html ). However, without implementation of the iout=2 option, these examples cannot be replicated using scipy.  Some popular sources such as Press et al.'s "Numerical Recipes" and Boost's odeint have solvers based on DOPRI5/DOP853 and provide access to dense output. 

** Details **

The original DOPRI5/DOP853 codes provide the option for a user-specified callback function to be called after each internal integration step.  Currently this callback function is only passed the solution at the current step.  With the proposed change, the callback function will be provided (optionally) with the coefficients of an interpolating polynomial for the solution between the last and current step by using the "dense output" option (iout=2).

As far as the user of scipy.integrate is concerned, the main change is the addition of an optional "dense_components" keyword argument to the ".set_solout" method of classes "ode" and "complex_ode".  Behavior should be unchanged in the absence of this keyword.  In addition, I have put in some additional error handling code for the bug reported here:  
and made some minor changes (see git commit messages) that should not impact existing code.

The new "dense_components" keyword argument is used to specify the components for which dense output is desired.  If this keyword option is present, the callback function will now be passed the additional information required for interpolation.  A convenience function "dense_dop" has been provided in _ode.py for interpolation using this information.

The introduction of the "dense_dop" function seems a bit awkward --- suggestions for alternative approaches are welcome.  Most users will simply want to evaluate the interpolation function to determine the solution at points between the last and current step position.   My motivation for development of access to the dense output is an application in which I will use the coefficients of the interpolating polynomial directly, rather than the interpolated solution.  Access to the interpolating polynomial may also be useful for event location.

I've updated the relevant docstrings and added test code for this new option (I have only tested this code on my Linux system with Python 2.7).

Comments are welcome.  Thanks -- Jim.

[Shampine, 1997]: Lawrence F. Shampine and Mark W. Reichelt, The MATLAB ODE Suite, SIAM J. Sci. Comput., 18(1), 1–22. 1997 (22 pages) http://dx.doi.org/10.1137/S1064827594276424

[Shampine, 2003], L. F. Shampine, I. Gladwell, S. Thompson, "Solving ODEs with Matlab", 2003, Cambridge University Press.

[Hairer, 1993] Hairer et al., "Solving Ordinary Differential Equations, Nonstiff Problems",  Second Revised Edition, Springer, 1993.

SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
Matthew Brett | 15 Apr 00:12 2015

Mass 1D interpolation, advice


In brain imaging (my world) - we need to do 1D interpolation of large
datasets, specifically, we have data in the order of (X, Y) = (200,
4000), where we are interpolating at the same ~200 X locations for
each of the 4000 rows.

We need to do this with extrapolation at the boundaries, so I think we
can't use scipy.interpolate.interp1d which is otherwise well suited
for the task.

scipy.interpolate.InterpolatedUnivariateSpline does deal with
extrapolation, but only works with vectors.

I think, in order to make speed reasonable, I will need to loop over
the rows in compiled code.

My questions are:

a) how easy would it be to try and link some external Cython code
against the fitpack code in scipy?
b) if that is not easy, what can I do that will also be of most use to
other scipy users?


Ralf Gommers | 6 Apr 00:31 2015

Growth and sustainability of Scipy development


Here's a long story about something I've been thinking about for a while. Some individual parts have been discussed with a few people, but the complete story is only how I look at things at the moment. I make some concrete proposals that could be agreed on in this thread if they turn out to be uncontroversial; we can start separate threads as needed though.


*Problem Statement*
The growth of core developer and project admin capacity, while pretty strong over the past few years, is not quite keeping up with the growth in contributions we're seeing.

Identify this as an important issue and start addressing it more actively.

*Main Proposals*
- recruitment: a "become a core dev" mentoring program
- docs:
    - write "how to become a core developer" documentation
    - improve the developer guide further
    - finally merge the Scipy 1.0 roadmap
- communication: set up a private core dev mailing list
- funding:
    - we got some donations (excellent, very useful)
    - have an agreed upon strategy (written down) of what to do with funds
    - then aim at getting more donations and/or getting grants
- governance: we have very little formalized, add something lightweight that covers most bases

Long Version

*Problem Statement*
The number of open PRs is very slowly increasing, and stands now at ~110. Some of those PRs are stalled, but many simply need reviewing. The number of active core developers who review and merge PRs is also increasing, but evidently not at the rate that is necessary for PRs to be reviewed quickly.  The sometimes long delays in getting a PR merged have likely also discouraged further contributions from some new contributors.

Note that we did have a lot of growth in active core developers over the past
few years - I remember Fernando Perez pointing out to me at EuroSciPy'11 that the
Scipy bus factor was exactly 2, now it's more like 8.

Therefore I'd like to look for new core developers in a more pro-active way.
That starts by very clearly saying:

  We need your help!

To be more precise, we need both domain specialists that help review PRs for
a single module and core developers that have a broad understanding of how
things are done in Scipy and that take responsibility for signing off on the
quality of PRs and are merging them independently.

A "core dev mentoring program".  This can go quite quickly for people that are experienced Python programmers and are familiar with the Scipy community; it may be as little as a few email conversations or Skype calls on how to get started.  On the other end of the scale, I'd also like to help people who are new to the Scipy community or even open source development get started - as long as they have the right level of coding/communication skills of course. Then the road to commit rights will of course be longer, but as long as the ambition to get there is present then it's worth the effort (I know I certainly wasn't very experienced in contributing to open source projects at the time I became Scipy release manager, and could definitely have used some help).

What developer docs are missing or need significant improvements? Let's add those (I volunteer for writing duty). We have a reasonable set (links at end of email), but it's quite scattered. And what we don't have is "how to become a core developer" documentation, let's add that. Plus finally merge the 1.0 roadmap: https://github.com/scipy/scipy/pull/2908.

I'd like to set up a new core dev closed mailing list. Not with the purpose to do less in public, but simply for organizing the communication that we already have off-list in a very ad-hoc manner now. Typical topics: deciding on commit rights, discussions on ranking GSoC students, potential funding sources and what to do with such funding. I'm not yet sure what to propose for people with access to that list. We currently have 20 people with commit rights, of which about half are actually active. Currently decisions are mostly made between the active devs.

I'm starting to see more opportunities to use funding in effective ways to sustain or accelerate Scipy development. Main ones are organizing developer meetups / sprints, and targeted funding to talented individuals (students mainly) and for work on particularly hairy bottlenecks (example: Windows test&build). We have a donation page on scipy.org and have received a number of donations, but have not yet actively solicited donations because of (a) no fiscal sponsorship agreement with NumFOCUS signed yet, and (b) while we can use the money now, I'm uncomfortable to spend large amounts without some written down strategy on how and why we spend funds.

Scipy activity is currently very much centered around code, and other project aspects get relatively little attention. Big email debates about governance are often not very productive (more high-bandwidth communication is necessary here), but it would be helpful to keep making steps in this area.  Documenting processes, how decisions are
made, who decides on things like commit rights, resolving API discussions,
etc. is important.  That goes hand in hand with the developer docs mentioned
above.  The HACKING and MAINTAINERS documents together with the Numpy
developer guide (which also applies to Scipy) do answer quite a few of those
questions, but they're not exactly comprehensive.

I think we can liberally borrow from other projects on this front, and draft something reasonable that covers at least basic aspects of how decisions are made and about things like licensing, (non-)ownership of trademarks and domain names. This will likely be much easier to reach consensus on than starting from scratch. 
Request: I'd like to keep detailed discussion on the content of such a governance doc draft for later, and here just discuss/decide whether or not we want to take this on now and who wants to be involved in putting the first draft together.

SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
Ralf Gommers | 4 Apr 10:54 2015

0.16.0 release

Hi all,

It's about time to start preparing for the 0.16.0 release. Unusually, there are zero blocking issues it looks like. A lot of PRs to merge of course, but nothing very large that still needs a lot of work.

Most important question: does anyone want to take the release manager role for this release? It would be good to start rotating this role for every minor release. This will help spread the knowledge on how to do releases and also spread the workload (it's not that much for a single release, but it adds up).


SciPy-Dev mailing list
SciPy-Dev <at> scipy.org

SciPy 2015 Conference Updates - LAST CALL FOR TALKS - extended to 4/10, registration open, keynotes announced, John Hunter Plotting Contest




SciPy 2015 will include 3 major topic tracks and 7 mini-symposia tracks. Submit a proposal on the SciPy 2015 website: http://scipy2015.scipy.org. If you have any questions or comments, feel free to contact us at: scipy-organizers <at> scipy.org. You can also follow <at> scipyconf on Twitter or sign up for the mailing list on the website for the latest updates!


Major topic tracks include:

- Scientific Computing in Python (General track)

- Python in Data Science

- Quantitative Finance and Computational Social Sciences


Mini-symposia will include the applications of Python in:

- Astronomy and astrophysics

- Computational life and medical sciences

- Engineering

- Geographic information systems (GIS)

- Geophysics

- Oceanography and meteorology

- Visualization, vision and imaging




Please register ASAP to help us get a good headcount and open the conference to as many people as we can. PLUS, everyone who registers before May 15 will not only get early bird discounts, but will also be entered in a drawing for a free registration (via refund or extra)! Register on the website at http://scipy2015.scipy.org




Keynote speakers were just announced and include Wes McKinney, author of Pandas; Chris Wiggins, Chief Data Scientist for The New York Times; and Jake VanderPlas, director of research at the University of Washington's eScience Institute and core contributor to a number of scientific Python libraries including sci-kit learn and AstroML.




In memory of John Hunter, creator of matplotlib, we are pleased to announce the Third Annual SciPy John Hunter Excellence in Plotting Competition. This open competition aims to highlight the importance of quality plotting to scientific progress and showcase the capabilities of the current generation of plotting software. Participants are invited to submit scientific plots to be judged by a panel. The winning entries will be announced and displayed at the conference. John Hunter’s family is graciously sponsoring cash prizes up to $1,000 for the winners. We look forward to exciting submissions that push the boundaries of plotting! See details here: http://scipy2015.scipy.org/ehome/115969/276538/ Entries must be submitted by April 13, 2015 via e-mail to plotting-contest <at> scipy.org




--Sprint, Birds of a Feather, Financial Aid and Talk submissions are open NOW

--Apr 10, 2015: Talk and Poster submission deadline

--Apr 13, 2015: Plotting contest submissions due

--Apr 15, 2015: Financial aid application deadline

--Apr 17, 2015: Tutorial schedule announced

--May 1, 2015: General conference speakers & schedule announced

--May 15, 2015 (or 150 registrants): Early-bird registration ends

--Jun 1, 2015: BoF submission deadline

--Jul 6-7, 2015: SciPy 2015 Tutorials

--Jul 8-10, 2015: SciPy 2015 General Conference

--Jul 11-12, 2015: SciPy 2015 Sprints


SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
James Phillips | 31 Mar 18:29 2015

Github code for Differential Evolution reevaluations

I have just reviewed the Differential Evolution code in github.  It is is particularly well coded and most algorithmically thorough, my professional compliments.

I noticed in the older DE code I use there is a possible performance improvement, and this seems to be the case in the github code as well.  Not every population vector changes within a generation, yet my older code reevaluates the unchanged vectors.  Of course this does not change the final results, but is a point of potential savings in total computation time.


James Phillips

SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
Andreea Georgescu | 31 Mar 04:34 2015

Pull request #4648


I was wondering if there's anything else I could improve on for the pull request #4688 (Fixes #4408: Vector-valued constraints in minimize() et al). It's still marked as "needs-work", although I followed all the things that Pauli Virtanen suggested. I'm a new contributor to SciPy so I'm not sure what the process is at this stage :). Should I just wait for it to be reviewed and merged, or is there anything else I should do?


SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
Fukumu Tsutsumi | 26 Mar 16:16 2015

About Application for GSoC

Hello, my name is Fukumu Tsutsumi. I'm a student at the University of Tokyo, Japan.
I'm planning to apply for GSoC with one of SciPy projects, but I deeply regret that I did not notice the deadline is so close (only about a day).
I wonder if I can meet the deadline. I'm very eager to participate in GSoC, so wouldn't like to give up this opportunity.
I would appreciate if anyone reply for me.


Fukumu Tsutsumi
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
Ralf Gommers | 25 Mar 20:19 2015

Re: Regarding taking up project ideas and GSoC 2015

On Wed, Mar 25, 2015 at 6:59 PM, Maniteja Nandana <maniteja.modesty067 <at> gmail.com> wrote:
Hi everyone,

I wanted to get some feedback on the application format and whether the mentioning of methods, API and other packages is necessary in the application or would it be preferable to provide a link to the Wiki page which contains that information.

Your proposal is already a lot more detailed than other proposals are, so I suggest to at least not make it any longer. Moving/keeping some of the background content to/in the wiki and linking to it in your proposal would be even better.
I would also update the timeline as early as possible, after I refine the ideas. It would also be great to have any other feedback.

Your timeline is now still empty, it's important to fill that in asap. It's easier to comment on a draft and improve it than suggest something from scratch. There are a number of things that have been suggested and you could put in (and a few I just thought of):

- write set of univariate test functions with known first and higher order derivatives
- same exercise for multivariate test functions
- define desired broadcasting behavior and implement
- refactor numdifftools.core.Derivative
- finalize API in a document
- integrate module into Scipy
- replace usages of numpy.diff with new scipy.diff functionality within Scipy
- (bonus points for at the end):
    - write a tutorial section about scipy.diff
    - write a nice set of benchmarks

On Mon, Mar 23, 2015 at 6:42 AM, Maniteja Nandana <maniteja.modesty067 <at> gmail.com> wrote:
Hi everyone, 
I was thinking it would be nice to put forward my ideas regarding the implementation of the package.

Thanks to Per Brodtkorb for the feedback.

On Thu, Mar 19, 2015 at 7:29 PM, <Per.Brodtkorb <at> ffi.no> wrote:



For your information I have reimplemented the approx._fprime and approx._hess code found in statsmodels and added the epsilon extrapolation

method of Wynn. The result you can see here:



This is wonderful, The main aim now is to find a way to determine whether the function is analytic, which is the necessity for the complex step to work. Though differentiability is one of the main necessities for analyticity, it would be really great if any new suggestions are there ?

I have also compared the accuracy and runtimes for the different alternatives here:



Thanks for the information. This would help me better in understanding the pros and cons for various methods.


Personally I like the class interface better than the functional one because you can pass the resulting object as function to other methods/functions and these functions/methods do not need to know what it does behind the scenes or what options are used. This simple use case is exemplified here:


>>> g = lambda x: 1./x

>>> dg = Derivative(g, **options)

>>> my_plot(dg)

>>> my_plot(g)


In order to do this with a functional interface one could wrap it like this:


>>> dg2  = lambda x: fprime(g, x, **options)

>>> my_plot(dg2)


If you like the one-liner that the function gives, you could call the Derivate class like this


>>> Derivate(g, **options)(x)


Which is very similar to the functional way:

>>> fprime(g, x, **options)

This is a really sound example for using classes. I agree that classes are better than functions with multiple arguments, and also the Object would e reusable for other evaluations.


Another argument for having it as a class is that a function will be large and

large functions are where classes go to hide”. This  is a quote of Uncle Bob’s that we hear frequently in the third and fourth Clean Coders episodes. He states that when a function starts to get big it’s most likely doing too much— a function should do one thing only and do that one thing well. Those extra responsibilities that we try to cram into a long function (aka method) can be extracted out into separate classes or functions.


The implementation in https://github.com/pbrod/numdifftools/blob/master/numdifftools/nd_cstep.py is an attempt to do this.


For the use case where n>=1 and the Richardson/Romberg extrapolation method, I propose to factor this out in a separate class e.g. :

>>> class NDerivative(object):

….      def __init__(self, f, n=1, method=’central’, order=2, …**options):


It is very difficult to guarantee a certain accuracy for derivatives from finite differences. In order to get error-estimates for the derivatives one must do several functions evaluations. In my experience with numdifftools it is very difficult to know exactly which step-size is best. Setting it too large or too small are equally bad and difficult to know in advance. Usually there is a very limited window of useful step-sizes which can be used for extrapolating the evaluated differences to a better final result. The best step-size can often be found around (10*eps)**(1./s)*maximum(log1p(abs(x)), 0.1) where s depends on the method and derivative order.  Thus one cannot improve the results indefinitely by adding more terms. With finite differences you can hope the chosen sampling scheme gives you reasonable values and error-estimates, but many times, you just have to accept what you get.  


Regarding the proposed API I wonder how useful the input arguments epsabs, epsrel  will be?

I was just then tinkering about controlling the absolute and relative errors of the derivative, but now it seems like we should just let the methods to take care of it.

I also wonder how one can compute the outputs abserr_round, abserr_truncate accurately?

This idea was from the implementation in this function. I am not sure of how accurate the errors would be, but I suppose this is possible to implement.



Best regards

Per A. Brodtkorb

Regarding the API, after some discussion, the class implementation would be something like


     Def __init__(f, h=None, method=’central’, full_output=False)

     Def __call__(self, x, *args, **kwds)



     Def __init__(f, h=None, method=’central’, full_output=False)

     Def __call__(self, x, *args, **kwds)



    Def __init__(f, h=None, method=’central’, full_output=False)

     Def __call__(self, x, *args, **kwds)



     Def __init__(f, h=None, method=’central’, full_output=False)

     Def __call__(self, x, *args, **kwds)


    Def __init__(f, n=1, h=None, method=’central’, full_output=False, **options)

    Def __call__(self, x, *args, **kwds)


Where options could be

Options = dict(order=2, Romberg_terms=2)

I would like to hear opinion on this implementation, where the main issues are

  1. whether the h=None default would mean best step-size found using by around (10*eps)**(1./s)*maximum(log1p(abs(x)), 0.1) where s depends on the method and derivative order or StepGenerator, based on  epsilon algorithm by wynn.
  2. Whether the *args and **kwds should be in __init__ or __call__, the preference by Perk was for it being in __call__ makes these object compatible with scipy.optimize.minimize(funx0args=()method=Nonejac=Nonehess=None,…..) where the args are passed both to the function and jac/hess if they are supplied.
  3. Are the input arguments for the __init__ sufficient ?
  4. What should we compute and return for full_output=True, I was thinking of the following options :

x: ndarray solution array,

success : bool a flag indicating if the derivative was calculated successfully        message : str which describes the cause of the error, if occurred            nfev int number of function evaluations 

abserr_round  : float absolute value of the roundoff error, if applicable            abserr_truncate : float absolute value of the truncation error, if applicable

It would be great any other opinions and suggestions on this.


SciPy-Dev mailing list
SciPy-Dev <at> scipy.org




SciPy-Dev mailing list
SciPy-Dev <at> scipy.org

SciPy-Dev mailing list
SciPy-Dev <at> scipy.org

SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
James Phillips | 24 Mar 20:00 2015

Robert Kern's accursed excellent DE implementation

I have been using Robert Kern's implementation of the Differential Evolution (DE) genetic algorithm for a decade, for the purpose of guessimating initial parameter estimates for curve fitting and surface fitting in my oprn source pyeq2 fitting library.

I can't find anything that works better, which gives rise to my current problem.

In trying to improve performance of my fitting library, I tried to use GPU calculations for each generation of the genetic algorithm, I found the following:

1) Robert's 2005 implementation of DE is not parallelizeable, as each crossover withing a generation can affect the population from which new items will be created.  That is, within a given generation the population *changes* as the algorithm runs the generation itself, and it must run serially in it's present form.

2) I can rework the algorithm to be parallelizeable by separating out "crossover", "breeding" and "evolving" into three separate steps, but months of testing show that population size and number of generations must beconsiderably  increased to match the results from Robert's version.  That is, making the algorithm parallelizable means slowing it down so I can speed it up!

I would like to increase performance, but cannot find any way to equal Robert's results without reducing performance prior to parallelization.

Any suggestions?

SciPy-Dev mailing list
SciPy-Dev <at> scipy.org
Ian Henriksen | 24 Mar 01:22 2015

GSOC 2015 projects

Hi all,

I'm putting together an application for GSOC 2015, and, although the deadline is approaching fast, I would still appreciate your feedback on a few of my ideas.

I am a masters student studying mathematics at Brigham Young University. Thus far my primary contribution to SciPy has been the Cython API for BLAS and LAPACK (see https://github.com/scipy/scipy/pull/4021).

My research is in Isogeometric Analysis, i.e. finite element analysis on spline curves. My open source work and research at BYU have given me a great deal of experience with Python and Cython, but I am also familiar with C, C++, and Fortran. As such, I have been reflecting on projects that would be best suited to my skill set, as well as most beneficial to SciPy.

I'm curious to know which of the following projects would be of greatest interest to the SciPy community:

1. Wrapping the PROPACK library for sparse linear algebra and using it for sparse SVD computation in scipy.sparse. There has been some initial work on f2py wrappers for PROPACK at https://github.com/jakevdp/pypropack, though it appears the wrappers are still incomplete.
2. Implementing an improved library for spline operations in SciPy. I am very familiar with the different refinement techniques used in CAD (knot insertion, degree elevation, etc.) and could implement a library that would be able to perform them all. My ideal here would be to write a C++ or Fortran library to do this and then wrap it via Cython. The emphasis would be primarily on writing code for refinement and evaluation that is both fast and general. I could include code for spline subdivision methods as well.
3. Adding support for Cython to both f2py and numpy.distutils. The goal here would be to allow f2py to generate cython-compatible wrappers from existing pyf files. I would also modify numpy.distutils so it could compile Cython files.
4. Wrap ffts (https://github.com/anthonix/ffts) and use it as an alternative to FITPACK in scipy.fft for use cases where it is faster.

Which of these projects would be most appreciated? I certainly want to be able to make a valid and, more importantly, useful contribution.


- Ian Henriksen
SciPy-Dev mailing list
SciPy-Dev <at> scipy.org