Mailman 3 February 2018 - SciPy-Dev

SciPy 1.0 paper writing proposal
by Ralf Gommers Jan. 29, 2019

Jan. 29, 2019

Hi all, TL;DR, let's write the long journal paper on SciPy that we've wanted for a while, let's form a small committee to coordinate, and get it out the door in 2-3 months. Motivation --------------- (credits for most of this text: Evgeni) Many scipy contributors' day jobs are in academia. Bibliometry -- papers in refereed journals and citations of papers by other papers -- is one of the main performance indicators in most academic establishments. Since we do not generate papers, scipy contributions are all but invisible for the purposes of a contributor's annual report. Of course, details vary wildly; in many cases a contributor manages to balance their time, or to argue common sense with their superiors, or get an approval for scipy work, or just ignores the issue altogether -- but sooner or later there is a form to be filled or boxes to be checked, and scipy contributions simply do not fit in. A peer-reviewed journal paper on scipy will help contributors get the academic credit they deserve. We can write *the* paper for SciPy 1.0, with overall project structure, goals, etc., and for specific features/modules a focus on say the last 3 years. History ---------- For SciPy 1.0 we had three targets on the publicity/credits front: an interesting release announcement, interesting blogs/stories (NumFOCUS blog, Hacker News, etc.) and a paper. We didn't have the bandwidth for a paper in the end, the rest was successful. [1] is a previous announcement on this list about writing (a) paper(s) on SciPy. We wanted both "short papers" to cover one or two releases (target journal JOSS) and a full paper as the authoritative reference for SciPy. We had an earlier attempt for a "short paper", it's mostly written but has stalled (see [2]). We ran out of steam on that one. To avoid that this time around, it would be good to have a clear public plan, target dates, and a small committee rather than one person to drive things forward. Proposal ------------ Here's a proposal for all aspects of this exercise that I can think about now. Some parts stolen from the AstroPy paper [3] (because their process worked quite well). Form a small coordination committee of 3-5 people that set up the paper structure, move things along when parts stall, propose/take decisions as needed, invite co-authors, and organise paper submission/rework. Paper writing to be done by whoever volunteers for a section, not just the coordination committee. First outline/structure to be done by committee, which then asks for review of structure and volunteers for section writing. Scope: a 6-10 page paper, covering history, package scope and structure, community/organisational aspects, key features and recent enhancements per module, and roadmap. Authorship: anyone who made a substantial contribution in the history of the project. Here "substantial" is interpreted as anything beyond a one-line doc fix. Rationale: better to be too inclusive than exclusive. Sign-up via a web form, we send the link to that form to all email addresses in the commit history till v1.0. Author order (details tbd by committee): 1. The SciPy Developers 2. Maintainers, paper writers, other key contributors - in order of contribution level 3. All other authors - alphabetically ordered Submission target: mid-April, to either PeerJ Computer Science or Journal of Open Research Software (tbd by committee). Comments? Volunteers for committee? References ---------------- [1] https://mail.python.org/pipermail/scipy-dev/2016-August/021474.html [2] https://github.com/scipy/scipy-articles/pull/4 [3] https://github.com/astropy/astropy-v2.0-paper Cheers, Ralf

12 19

General discussion on parallelisation
by Andrew Nelson Sept. 8, 2018

Sept. 8, 2018

As part of https://github.com/scipy/scipy/pull/8259 I'm proposing that a `workers` keyword is added to optimize.differential_evolution to parallelise some computation. The proposal is that: 1. the workers keyword accepts either an integer or an object with a map-like method. 2. If an integer is supplied then the parallelisation is taken care of by scipy (more on that later), with -1 signifying that all processors are to be used. 3. If an object with a map-like method is supplied, e.g. `multiprocessing.Pool.map`, `mpi4py.futures.MPIPoolExecutor.map`, etc, then the parallelisation is taken care of by that object. This allows the user to specify the parallelisation configuration for their problem. 4. If workers=1, then computation will be done by the builtin `map` function. Now we come to the under the hood part. I've written something called PoolWrapper ( https://github.com/andyfaff/scipy/blob/b14bb513c0ffb9807a67663d39b9ab399375…) which wraps `multiprocessing.Pool` to achieve the behaviour outlined above. It can be used as a context manager, or the user of the object can decide when to close the resources opened by PoolWrapper. I've looked at using joblib instead of PoolWrapper and it seems useful but it doesn't have a couple of bits of functionality that are needed for this specific problem: viz: 5. joblib.Parallel doesn't have a map method (desirable to allow 3) so a small wrapper would have to be created anyway. 6. joblib.Parallel creates/destroys a multiprocessing.Pool each time the Parallel object is `__call__`ed. This leads to significant overhead. One can use the Parallel object with a context manager, which allows reuse of the Pool, but I don't think that's do-able in the context of using the DifferentialEvolutionSolver (DES) object as an iterator: >>> solver = DifferentialEvolutionSolver(func, bounds) >>> # use DES object as an iterator >>> for it in solver: ... res = next(solver) >>> print(res) >>> # use the DES.solve method >>> res = solver.solve() Whilst the DES object is not currently public (it's called by the differential_evolution function) it would be nice to expose it in the future, and people will want to use both approaches. Unfortunately with the first approach if we used joblib.Parallel we'd have to use Parallel.__call__ in DES.next() which has the overhead of creating/destroying Pools. For efficient use of resources the Pool should persist for the lifetime of the DES object. I also looked into `concurrent.futures.ProcessPoolExecutor`, but it's not available for Python 2.7. The purpose of this email is to elicit feedback for developing parallelisation strategy for scipy - what does the public interface look like, what does scipy do under the hood? Under the hood I think a mixture of PoolWrapper and joblib.Parallel could be used (with scipy vendoring joblib). A. -- _____________________________________ Dr. Andrew Nelson _____________________________________

5 18

Are you using Bento to build scipy?
by CJ Carey July 16, 2018

July 16, 2018

With the 1.0 release coming up, I figured this is a good time to take the pulse of the community regarding Bento. The project seems to have been abandoned before reaching stability, and it seems that the core problems it was trying to address have been at least partially fixed by the more official build systems. If you use Bento and would like to keep using it post 1.0, or if my previous sentence is wrong, please let me know!

10 15

Stochastic calculus package
by Maurizio Cipollina June 15, 2018

June 15, 2018

Hi everybody, and first of all thanks for the great job you have been doing, on which I relied *A LOT* over time. Stochastic calculus and Stochastic Differential Equations are a significant branch of mathematics that to my understanding is underrepresented in scipy, as well as in the open source python stack at large. I may well be missing something, but as I can see one can find out there many recipies or collections of recipies, and a number of packages/projects skewed towards finance, but no structured framework for generating and handling stochastic processes and numerically solving SDEs. This is a pity, since scipy/numpy provide efficen implementations of all the basic ingredients needed to do the trick. Out of professional need, and prompted by this gap, I have developed a general purpose package of which I own the copyright, and which I would be happy to release under the BSD license as part of scipy, in case there is a consensus among the scipy community that it makes sense to do so (the package has no dependencies beyond standard library, numpy and scipy, and might probably fit best as a subpackage of scipy.stats). Before setting up a PR, I would be happy to hear your thoughts: I have pasted below the package docstring, that should give an overview of its scope and limitations. Please note that some of the realized processes stray into mathematical finance territory, somehow inevitably since finance is one major field of application of stochastic calculus, but the focus is the latter, not the former, at least in my intentions. Thanks in advance for taking your time to go through this, and for your comments and suggestions. Maurizio “”” =========================================================== Stochastic - Monte Carlo simulation of stochastic processes =========================================================== This package provides: 1. A `process` class, a subclass of `numpy.ndarray` representing a sequence of values in time, realized in one or several paths. Algebraic manipulations and ufunc computations are supported for instances that share the same timeline and comply with numpy broadcasting rules. Interpolation along the timeline is supported via callability of `process` instances. Process-specific functionalities, such as averaging and indexing along time or across paths, are delegated to process-specific methods, attributes and properties (no overriding of `numpy.ndarray` operations). 2. The `source` class and its subclasses, stochasticity sources providing numerical realizations of the differentials commonly found in stochastic differential equations (SDEs), with or without memory of formerly invoked realizations. 3. The `integrator` class, a general framework for numerical SDE integration, intended for subclassing, to assist the definition of step by step integration algorithms and the computation of numerical realizations of stochastic processes. 4. Several objects providing realizations of specific stochastic processes along a given timeline and across a requested number of paths. Realizations are obtained via explicit formulae, in case the relevant SDE has a closed form solution, or otherwise as `integrator` subclasses performing Euler-Maruyama numerical integration. 5. A `montecarlo` class, as an aid to cumulate the results of several Monte Carlo simulations of a given stochastic variable, and to extract summary estimates for its probability distribution function and statistics, thus supporting simulations across a number of paths that exceeds the maximum allowed by available memory. 6. Several analytical results relating to the supported stochastic processes, as a general reference and for testing purpouses. For all sources and realizations, process values can take any shape, scalar or multidimensional. Correlated multivariate stochasticity sources are supported. Time-varying process parameters (correlations, intensity of Poisson processes, volatilities etc.) are allowed whenever applicable. `process` instances act as valid stochasticity source realizations (as does any callable object complying with a `source` protocol), and may be passed as a source specification when computing the realization of a given process. Computations are fully vectorized across paths, providing an efficient infrastructure for simulating a large number of process realizations. Less so, for large number of time steps: integrating 1000 time steps across 100000 paths takes seconds, one million time steps across 100 paths takes minutes. General infrastructure ====================== .. autosummary:: :toctree: generated/ process -- Array representation of a process (a subclass of -- numpy.ndarray). kfunc -- Base class for functions with managed keyword arguments. source -- Base class for stochasticity sources. true_source -- Base class for stochasticity sources with memory. integrator -- General framework for numerical SDE integration. montecarlo -- Summmary statistics of Monte Carlo simulations. Stochasticity sources ===================== .. autosummary:: :toctree: generated/ wiener -- dW, a source of standard Wiener process (Brownian -- motion) increments. symmetric_wiener -- as above, with symmetric paths (averages exactly 0). true_wiener -- dW, a source of standard Wiener process (Brownian -- motion) increments, with memory. poisson -- dN, a source of Poisson process increments. compound_poisson -- dJ, a source of compound Poisson process increments. -- (jumps distributed in time as a Poisson process, -- and in size as a given `scipy.stats` distribution). compound_poisson_rv -- Preset jump size distributions for compound Poisson -- process increments. Stochastic process realizations =============================== .. autosummary:: :toctree: generated/ const_wiener_process -- Wiener process (Brownian motion), with -- time-independent parameters. const_lognorm_process -- Lognormal process, with time-independent -- parameters. wiener_process -- Wiener process (Brownian motion). lognorm_process -- Lognormal process. ornstein_uhlenbeck_process -- Ornstein-Uhlenbeck process (mean-reverting -- Brownian motion). hull_white_process -- F-factors Hull-White process (sum of F -- correlated mean-reverting Brownian -- motions). hull_white_1factor_process -- 1-factor Hull-White process (F=1 Hull-White -- process with F-index collapsed to a -- scalar). cox_ingersoll_ross_process -- Cox-Ingersoll-Ross mean-reverting process. full_heston_process -- Heston stochastic volatility process -- (returns both process and volatility). heston_process -- Heston stochastic volatility process -- (stores and returns process only). jump_diffusion_process -- Jump-diffusion process (lognormal process -- with compound Poisson jumps). merton_jump_diffusion_process -- Merton jump-diffusion process -- (jump-diffusion process with normal -- jump size distribution). kou_jump_diffusion_process -- Kou jump-diffusion process (jump-diffusion -- process with double exponential -- jump size distribution). Shortcuts:: lognorm -- lognorm_process oruh -- ornstein_uhlenbeck_process hwp -- hull_white_process hw1f -- hull_white_1factor_process cir -- cox_ingersoll_ross_process heston -- heston_process mjd -- merton_jump_diffusion_process kou -- kou_jump_diffusion_process Analytical results ================== Exact statistics for the realized stochastic processes are listed below, limited to the case of constant process parameters and with some further limitations to the parameters' domains. Function arguments are consistent with those of the corresponding processes. Suffixes `pdf`, `cdf` and `chf` stand respectively for probability distribution function, cumulative probability distribution function, and characteristic function. Black-Scholes formulae for the valuation of call and put options have been included (prefixed with `bs` below). .. autosummary:: :toctree: generated/ wiener_mean wiener_var wiener_std wiener_pdf wiener_cdf wiener_chf lognorm_mean lognorm_var lognorm_std lognorm_pdf lognorm_cdf lognorm_log_chf oruh_mean oruh_var oruh_std oruh_pdf oruh_cdf hw2f_mean hw2f_var hw2f_std hw2f_pdf hw2f_cdf cir_mean cir_var cir_std cir_pdf heston_log_mean heston_log_var heston_log_std heston_log_pdf heston_log_chf mjd_log_pdf mjd_log_chf kou_mean kou_log_pdf kou_log_chf bsd1d2 bscall bscall_delta bsput bsput_delta Notes ===== To the benefit of interactive and notebook sessions, and as an aid to fluently freeze or vary the plurality of parameters that define each stochastic process, all sources, process realizations and analytical results are objects with managed keywords (subclasses of `kfunc`): if the corresponding classes are instantiated with both variables (non keyword arguments) and parameters (keyword arguments), they behave as functions returning the computations' result; if called with parameters only, return an instance that stores the set parameters, and exposes the same calling pattern (call with variables and optionally with parameters to get results, call with parameters only to get a new instance storing modified parameters). Parameters that are not specified fall back to the class defaults, if calling the class, or to the instance's stored values, if calling an instance. """

4 7

Using Pythran to compile some of the scipy internals
by Serge Guelton May 11, 2018

May 11, 2018

Hi Scipy-dev, = Context I've started the experiment of compiling some of the scipy kernels with the Pythran[0] compiler. The approach was basically to pick some easy .pyx file, turn them back to high-level Python + Numpy files then processing them with Pythran. This results in easier to maintain, Python-compatible files that run at the same speed as their Cython alternative up to 10 times faster depending on the kernels and whether SIMD instructions are activated or not (I did not play with parallelization in that experiment). This required a few improvement to Pythran, but not that much, see https://github.com/serge-sans-paille/pythran/pull/770 for the associated PR on pythran, which includes the following scipy kernels: hausdorff.py max_len_seq_inner.py solve_toeplitz.py spectral.py = Why bother? According to that first round of experiments, using Pythran can bring extra speed improvements while remaining at a higher level than Cython. Both points looks like an improvement to me. = Cost? This adds an extra dependency on Pythran, which uses C++ as backend. This increases the failure surface. Although alive since 2012 and being tested a lot [1] on Linux (but scarcely on Windows), its is obviously less mature than cython = Alternatives There is an experimental Pythran mode in cython[2] that uses Pythran as a backend for numpy operations. Unfortunately it is still at early stages and cannot translate calls to ``np.roll`` or ``np.sum(a[i,:] - b[j, :])`` while Pythran supports it. Instead of translating Cython files, I could also focus on some pure-python functions. I tested Pythran on the rosenbrock function and I get good speedup (from 1.5x to 4x depending on vectorization being enabled or not) there too. So yeah, that's a rather long introduction to probe the interest here around that idea :-) Best, Serge PS: I started https://github.com/scipy/scipy/pull/8306 meanwhile, but let's discuss here first :-) [0] https://github.com/serge-sans-paille/pythran [1] https://travis-ci.org/serge-sans-paille/pythran/builds/330052310 [2] http://cython.readthedocs.io/en/latest/src/userguide/numpy_pythran.html

8 14

GSoC 2018: Cythonizing
by Lars G. March 7, 2018

March 7, 2018

Dear SciPy devs, I'm currently thinking about an application for this year's GSoC as well. As there already seems to be a large interest in the rotation formalism I'm trying to find another area that matches my interest and skill. I've dug up this proposal in scikit-image from GSoC 2015 https://github.com/scikit-image/scikit-image/wiki/GSoC-2015#rewriting-scipy… and judging by the state of scipy/ndimage/src/ nobody has worked on this proposal yet (feel free to correct me). Alternatively I could imagine something similar for other sub-packages, e.g. scipy/signal which features many source files in C as well. So basically if there is an interest I could try to port C / Python code to Cython. What I would like to know: - Is there an interest? ;) - Is the original proposal in scikit-image still unfinished and are the potential mentors still interested in mentoring? - If there is a general interest to cythonize C or Python code during a GSoC project, which parts / sub-packages of SciPy would you priorize? As for my current involvement with SciPy: - I've already added a small function written in Cython https://github.com/scipy/scipy/pull/8350 - as part of a larger PR extending the signal module https://github.com/scipy/scipy/pull/8264 which will possibly merged this week. - I already cythonized slow parts of the above PR and plan to add these with new PRs after #8264 is merged. If this receives positive feedback I'd be happy to draft a more complete proposal / application based on the discussion around this. Best regards, Lars

7 12

WIP: Class based Optimizers
by Andrew Nelson March 6, 2018

March 6, 2018

The #8414 PR on github (https://github.com/scipy/scipy/pull/8414) introduces an `Optimize` class for `scipy.optimize`. The `Optimize` class can be advanced like an iterator, or be advanced to convergence with the `solve` method. The PR is intended to provide a cleaner interface to the minimisation process, with the user being able to interact with the solver during iteration. This allows the variation of solver hyper-parameters, or other user definable interactions. As such `callback`s become redundant, as you can access all the optimizer information required. Other historical parameters may become redundant in these objects, such as `disp`. This is also a chance for a ground-zero approach to how people interact with optimizers. The PR is intended as a Request For Comment to gain feedback on architecture of the solvers. Given that it's a key part of infrastructure this may create a large amount of feedback, if I don't reply to it all bear in mind that I don't have unlimited time resource. In it's current state I've implemented `Optimizer`, `Function`, `LBFGS` and `NelderMead` classes. `LBFGS` and `NelderMead` are being called by the corresponding functions already in existence. The current test suite passes, but no unit tests have been written for the functionality. I'd prefer to add such a specific test suite once the overall architecture is settled. These optimisers were straightforward to translate. However, I don't know how this approach would translate to other minimisers, or things like `least_squares`. If a design could be agreed, then perhaps the best approach would be to tackle the remaining optimizers on a case by case basis once this PR is merged.

9 16

Spline interpolation in ndimage
by Jaime Fernández del Río March 6, 2018

March 6, 2018

Hi all, We have been discussing in issue #8465 <https://github.com/scipy/scipy/issues/8465> about the extension modes for the interpolation module of ndimage. As some other issues linked in that one show, this is an old problem, that has been recurringly discussed, and there mostly is agreement on what the correct behavior should be. But as all things ndimage, implementing any change is hard, in this case not only because the code is complicated, but also because of the mathematical subtleties of the interpolation method used. In trying to come up with a way forward for this I think I can handle the code complexity, but am having trouble being sure that I come up with a sound math approach. I think I have more or less figured it out, but I don't have a good academic background on signal processing, so was hoping that someone who does (I'm thinking of you, Chuck!) can read my longish description of things below and validate or correct it. My main source of knowledge has been: M. Unser, "Splines: a perfect fit for signal and image processing," in IEEE Signal Processing Magazine, vol. 16, no. 6, pp. 22-38, Nov 1999. Recommendations for bigger, better readings on the subject are also welcome! ----- In what follows I'll assume the input image is 2-D and NxN for simplicity, but I think everything generalizes easily If I'm understanding it right, all of our interpolation functions use B-splines for interpolation. So instead of having NxN discrete values on a grid, we have NxN B-splines centered at the grid points, and we keep the NxN coefficients multiplying each spline to reproduce the original image at the grid points exactly. If the spline order is < 2, the spline coefficients are the same as the image values. But if order >= 2 (and we use 3 by default), these have to be calculated in a non-trivial fashion. This can be done efficiently by applying a separable filter, which is implemented in ndimage.spline_filter. Because B-splines have compact support, when using splines of order n we only need to consider the B-splines on an (n + 1)x(n + 1) grid neighborhood of the point being interpolated. This is more or less straightforward, until you move close to the edges and all of a sudden some of the points in your grid neighborhood fall outside the original image grid. We have our extend mode, which controls how this points outside should be mapped to points inside. But here is where things start getting tricky... When the spline coefficients are computed (i.e. when ndimage.spline_filter is called), assumptions have to be made about boundary conditions. But ndimage.spline_filter does not take a mode parameter to control this! So what I think ndimage does is compute the spline coefficients assuming "mirror symmetric" boundary conditions, i.e.: a b c d c b | a b c d | c b a b c d So if our interpolated point is within the image boundaries, but some of the grid points in its (n + 1)x(n + 1) neighborhood fall outside the boundary, the code uses a mirror symmetric extension mode to fill in those values. This approach has a funny smell to it, but even if it's formally incorrect I think it would only be marginally so, as the values are probably going to be very close to the correct ones. The problem comes when the point being interpolated is outside the boundaries of the image. We cannot use mirror-symmetric spline coefficients to extend if e.g. we have been asked to extend using wrap mode. So what ndimage does is first map the value outside the boundaries to a value within the boundaries, using the given extension mode, then interpolate it as before, using mirror-symmetric coefficients if needed because its (n + 1)x(n + 1) neighborhood extends outside. Again, this smells funny, but it is either correct or very close to correct. This is mostly all good and well, except for the "wrap" and "reflect" extension modes: in these cases the area within one pixel of the image boundaries is different from anything inside the image, so we cannot use that approach. So what ndimage does is basically make shit up and use something similar, but not quite right. "reflect" is mostly correct, except for within that pixel of the boundary, but "wrap" is a surprising and confusing mess. So how do we fix this? I see two ways forward: 1. What seems the more correct approach would be to compute the spline coefficients taking into account the extension mode to be used, then use the same extension mode to fill in the neighborhood values when interpolating for a point outside the boundaries. 1. First question is whether this is doable? I need to work out the math, but for "wrap" it looks like it should be, not 100% sure if also is for "reflect". 2. Assuming it is it has the main advantage of doing things in a more general and understandable way once you have enough background knowledge. 3. It does go a little bit against our API design: you can control whether the input is spline-filtered automatically with a parameter, the idea being that you may want to do the filtering yourself if you are going to apply several different transformations to the same image. If the mode of the filtering has to be synced with the mode of the transformation, letting the user do it themselves is a recipe for disaster, because it's going to lead to very hard to track bugs. 4. As elsewhere in ndimage, the current implementation does a lot of caching, which works because it always interpolates for a point within the image boundaries. If we started interpolating for points outside the boundaries without first mapping to within there may be a performance hit which has to be evaluated. 2. The other approach is, for "wrap" and "reflect" modes, pad the input image with an extra pixel in each direction, then compute our current "mirror symmetric" spline coefficients, and leave things as they are right now, aside from some changes to the mapping of values to take the extra pixels into account. 1. This looks like a nightmare of special cases everywhere and potential off-by-one errors while putting it together, but it would just go along with the ugliness of the existing code. 2. It's unclear what we would do if we are given an input with prefilter=False, so this also breaks the current API, probably even more so. Any thoughts or recommendations are very welcome! Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

2 3

Contributing to SciPy through GSoC
by Anubhav Patel March 4, 2018

March 4, 2018

Hi everyone, I want to work on SciPy as part of GSoC and I have few queries. 1. On the Ideas Page, there was a mention of scipy.spatial.transform module. I want to know what will be the exact purpose of this module? 2. Whether the idea for a module for numerical differentiation was dropped completely? 3. Apart from those ideas listed on ideas page, are there any other area where you guys would like to see contribution on?

4 4

SciPy IRC channel lacks moderation
by Mike Graham March 3, 2018

March 3, 2018

Hi, The SciPy IRC channel on Freenode has been without moderation for a while. On rare occasion, we get trolls in the channel. For example, recently we had a couple instances of a troll posting the n-word in channel and tagging many users (most likely with the intent of sending them notifications). It would also be nice if someone could do things like fix dead links in the channel's topic statement. I would suggest that either the channel be discontinued (and removed from https://www.scipy.org/scipylib/mailing-lists.html ) or that the SciPy organization appoint a "Group Contact" with Freenode to administrate the channel. (Freenode requires the organization a channel is named after to notify them of the desired GC.) Nathan Goldbaum and Nathaniel J. Smith are channel regulars -- if you are looking for someone to appoint as GC, you might try to recruit one of them. Best, Mike

5 9