Mailman 3 October 2021 - SciPy-Dev

creating a new scipy/archive repo
by Ralf Gommers Oct. 10, 2021

Oct. 10, 2021

Hi all, For things like logos, presentations, meeting minutes, and other SciPy project content that doesn't belong in the main repo, I'd like to make a new `archive` repo. This is analogous to https://github.com/numpy/archive, which has worked pretty well. The need for this came up again when looking for the SciPy logo, which is hard to find (there's links in some old issues). Once we have such a repo, I'm sure we'll find more things to put in it. I plan to make this repo on the weekend, please speak up if there are any concerns or better ideas. Cheers, Ralf

5 5

Re: Shifting `NumericalInverseHermite` and `rvs_ratio_uniforms` to UNU.RAN
by Christoph Baumgarten Oct. 7, 2021

Oct. 7, 2021

Hi all, just a few more comments on this topic: - The addition of NaiveRatioUniforms (PR 14717, [2]) does not lead to any issues. For now, we just keep the existing function rvs_ratio_uniforms and the user can choose which one to use. - Replacing the Python implementation of NumericalInversionHermite with the one from UNU.RAN was already planned when Matt initially implemented the functionality. Overall, the change happens in the background, i.e. the API is almost unchanged with the exception of the two points that Tirth mentioned: a) midpoint_error is set to nan and the deprecation message informs the user that there is a "better" error measure (the current value is just an approximation of the error that can be estimated in the new version) => the main question is whether this approach to deprecate midpoint_error is correct b) the deprecation of tol in favour of the new parameter u_resolution is generally backwards-compatible: tol can still be provided and will be preferred over u_resolution in that case. However, if neither of the values is specified, the new default is 1e-10 (to be consistent with NumericalInversePolynomial) instead of 1e-12. From a practical point of view, this should not make a difference: both errors are very small and differences in the generated rvs will be minimal. One could also keep the default at 1e-10, but even in that case one cannot guarantee that the values are exactly the same when switching from Python to C. Finally, there can be differences in the performance. However, if used properly (i.e., by avoiding to use rv_continuous objects), the setup times should be comparable overall while the actual sampling should be quite a bit faster with the new approach (that might be something we need to test a bit more). So Tirth and I hope that the approach taken in the PR is acceptable. Thanks for your feedback. Christoph On Wed, Oct 6, 2021 at 6:01 PM <scipy-dev-request(a)python.org> wrote: > Send SciPy-Dev mailing list submissions to > scipy-dev(a)python.org > > To subscribe or unsubscribe via the World Wide Web, visit > https://mail.python.org/mailman3/lists/scipy-dev.python.org/ > or, via email, send a message with subject or body 'help' to > scipy-dev-request(a)python.org > > You can reach the person managing the list at > scipy-dev-owner(a)python.org > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of SciPy-Dev digest..." > > Today's Topics: > > 1. Shifting `NumericalInverseHermite` and `rvs_ratio_uniforms` to > UNU.RAN > (Tirth Patel) > > > ---------------------------------------------------------------------- > > Message: 1 > Date: Wed, 6 Oct 2021 17:38:14 +0530 > From: Tirth Patel <tirthasheshpatel(a)gmail.com> > Subject: [SciPy-Dev] Shifting `NumericalInverseHermite` and > `rvs_ratio_uniforms` to UNU.RAN > To: scipy-dev <scipy-dev(a)python.org> > Message-ID: > <CABpuv395jjVZ6GrUCK17e=CUm5fUt= > WzMmEOLDp8FPsSZF2bqQ(a)mail.gmail.com> > Content-Type: multipart/alternative; > boundary="00000000000077a9ea05cdae02d5" > > Hi all, > > I and Christoph (@chrisb83) have been working on shifting some existing > random variate sampling methods to UNU.RAN. [1] shifts > `scipy.stats.NumericalInverseHermite` to a UNU.RAN based implementation and > [2] replicates `scipy.stats.rvs_ratio_uniforms` using the > `NaiveRatioUniforms` method in UNU.RAN. We have run into some breaking > changes/deprecations and would appreciate it if someone can help review/ > approve them. > > For `NumericalInverseHermite`, there are two big changes: > > 1. This method has a public attribute `midpoint_error` that reports the > approximation errors in the interpolating spline. But UNU.RAN doesn't > provide such an attribute or any functionality to somehow retrieve it. The > way we currently deal with this is to return `np.nan` and raise a > `DeprecationWarning` whenever `midpoint_error` is accessed. Instead of > suddenly removing `midpoint_error`, we decided to change its value and run > a deprecation cycle. > 2. The parameter `tol` has been deprecated in the favor of `u_resolution` > and the default value has been changed from `1e-12` to `1e-10`. > > One potential downside of moving to UNU.RAN is (very) slow setup with > `scipy.stats` distributions. > > On master: > > In [1]: from scipy.stats import norm, NumericalInverseHermite > > In [2]: dist = norm() > > In [3]: %timeit NumericalInverseHermite(dist) > 26.8 ms ± 1.03 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) > > With UNU.RAN: > > In [1]: from scipy.stats import norm, NumericalInverseHermite > > In [2]: dist = norm() > > In [3]: %timeit NumericalInverseHermite(dist, u_resolution=1e-12) > 2.46 s ± 35.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) > > One way to get over it is to avoid vectorizing the functions by using math > functions: > > In [1]: from math import pi, exp, sqrt, erf > > In [2]: class StandardNormal: > ...: def pdf(self, x): > ...: return 1/sqrt(2*pi) * exp(-0.5 * x*x) > ...: def cdf(self, x): > ...: return (1.0 + erf(x / sqrt(2.0))) / 2.0 > ...: > > In [3]: dist = StandardNormal() > > In [4]: %timeit NumericalInverseHermite(dist, u_resolution=1e-12) > 14.3 ms ± 922 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) > > As methods like this are primarily supposed to be used with custom > distributions, we can redirect the users to instead use the `rvs` method > for the distributions already implemented in `scipy.stats`. The tutorial > [3] clearly states this and the slowdown with `scipy.stats` distributions > is noted there. > > We also noted that sampling with UNU.RAN's `NaiveRatioUniforms` method is > about ~10x slower than `rvs_ratio_uniforms` because it calls a density > function implemented in Python which is (very) inefficient. This problem > can be addressed by adding support for `LowlevelCallable`s which I hope to > add in a follow-up PR. > > It would be very helpful to get any reviews/approval on the changes > mentioned above. Let us know what you think! > > [1]: https://github.com/scipy/scipy/pull/14703 > [2]: https://github.com/scipy/scipy/pull/14717 > [3]: https://scipy.github.io/devdocs/tutorial/stats/sampling.html > > > -- > Thanks, > Tirth Patel >

1 0

Krylov Subspace Method Scipy
by Fabian Schneider Oct. 7, 2021

Oct. 7, 2021

Hello, recently I implemented the Krylov Subspace Method for approximating the matrix exponential, respectively the phi-Functions, for very large sparse matrices in systems of ODEs. Paper by Saad for the normal exponential: https://www-users.cse.umn.edu/~saad/PDF/RIACS-90-ExpTh.pdf Paper by Niesen/Wright for phi-Functions: https://dl.acm.org/doi/10.1145/2168773.2168781 I would be happy to contribute my code, if you think that's a useful addition. Best regards Fabian Schneider

1 0

Shifting `NumericalInverseHermite` and `rvs_ratio_uniforms` to UNU.RAN
by Tirth Patel Oct. 6, 2021

Oct. 6, 2021

Hi all, I and Christoph (@chrisb83) have been working on shifting some existing random variate sampling methods to UNU.RAN. [1] shifts `scipy.stats.NumericalInverseHermite` to a UNU.RAN based implementation and [2] replicates `scipy.stats.rvs_ratio_uniforms` using the `NaiveRatioUniforms` method in UNU.RAN. We have run into some breaking changes/deprecations and would appreciate it if someone can help review/ approve them. For `NumericalInverseHermite`, there are two big changes: 1. This method has a public attribute `midpoint_error` that reports the approximation errors in the interpolating spline. But UNU.RAN doesn't provide such an attribute or any functionality to somehow retrieve it. The way we currently deal with this is to return `np.nan` and raise a `DeprecationWarning` whenever `midpoint_error` is accessed. Instead of suddenly removing `midpoint_error`, we decided to change its value and run a deprecation cycle. 2. The parameter `tol` has been deprecated in the favor of `u_resolution` and the default value has been changed from `1e-12` to `1e-10`. One potential downside of moving to UNU.RAN is (very) slow setup with `scipy.stats` distributions. On master: In [1]: from scipy.stats import norm, NumericalInverseHermite In [2]: dist = norm() In [3]: %timeit NumericalInverseHermite(dist) 26.8 ms ± 1.03 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) With UNU.RAN: In [1]: from scipy.stats import norm, NumericalInverseHermite In [2]: dist = norm() In [3]: %timeit NumericalInverseHermite(dist, u_resolution=1e-12) 2.46 s ± 35.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) One way to get over it is to avoid vectorizing the functions by using math functions: In [1]: from math import pi, exp, sqrt, erf In [2]: class StandardNormal: ...: def pdf(self, x): ...: return 1/sqrt(2*pi) * exp(-0.5 * x*x) ...: def cdf(self, x): ...: return (1.0 + erf(x / sqrt(2.0))) / 2.0 ...: In [3]: dist = StandardNormal() In [4]: %timeit NumericalInverseHermite(dist, u_resolution=1e-12) 14.3 ms ± 922 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) As methods like this are primarily supposed to be used with custom distributions, we can redirect the users to instead use the `rvs` method for the distributions already implemented in `scipy.stats`. The tutorial [3] clearly states this and the slowdown with `scipy.stats` distributions is noted there. We also noted that sampling with UNU.RAN's `NaiveRatioUniforms` method is about ~10x slower than `rvs_ratio_uniforms` because it calls a density function implemented in Python which is (very) inefficient. This problem can be addressed by adding support for `LowlevelCallable`s which I hope to add in a follow-up PR. It would be very helpful to get any reviews/approval on the changes mentioned above. Let us know what you think! [1]: https://github.com/scipy/scipy/pull/14703 [2]: https://github.com/scipy/scipy/pull/14717 [3]: https://scipy.github.io/devdocs/tutorial/stats/sampling.html -- Thanks, Tirth Patel

1 0

Proposal to add a binary feature coincidence test
by James Mathews Oct. 4, 2021

Oct. 4, 2021

I would like to propose a minor, stand-alone statistical test for inclusion in the stats submodule. It assesses the statistical significance of the number of subsamples of a given sample having all features out of a list of binary feature variables with given frequencies. I.e. "coincidence". The proposal seems to meet several of the expectations described in the developer guide for a feature addition: - "It fits the topic of the submodule, and does not require extensive support frameworks to operate." It is a statistical test. And the implementation is completely stand-alone. - "The implementation looks sound and unlikely to need much tweaking in the future." This will be up to reviewers to gauge, but the preliminary implementation, https://github.com/nadeemlab/coincidencetest/blob/main/coincidencetest/_coi… is documented, tested, https://github.com/nadeemlab/coincidencetest/blob/main/tests/test_coinciden… and largely passes the linter / code style checker. The underlying mathematics is written up here https://arxiv.org/abs/2109.13876 and was recently submitted for publication. - "Someone wants to contribute it." Check! However one of the indicated points admittedly requires some discussion: - "The method is applicable in many fields and 'generally agreed' to be useful." The test applies in the context of any number of binary variables, which is certainly applicable in many fields. However, because it is recent work, it is of course not yet generally agreed to be useful. It would be reasonable to object to inclusion on the grounds that the work is too recent and has not proven its utility in community usage, but we hope that this objection can be overcome by soliciting independent review. I'll try to answer the questions about this bullet point which are also indicated in the developer guide: - "Is the method used/useful in different domains in practice? How much domain-specific background knowledge is needed to use it properly?" We at the Nadeem lab (https://nadeemlab.org) have used the method in several bioinformatics-oriented situations. But it is just as applicable in engineering or social science contexts, and very little background is needed to apply it. - "Consider the code already in the module. Is what you are adding an omission? Does it solve a problem that you'd expect the module be able to solve? Does it supplement an existing feature in a significant way?" The tests portion of the stats submodule offers ways to quantify statistical significance in various situations. What we are adding is a high-order method (many variables), with a high degree of generality of application. It directly supplements the methods one might call second-order (in fact, by literally generalizing the Fisher exact test). - "Consider the equivalence class of similar methods / features usually expected. Among them, what would in principle be the minimal set so that there's not a glaring omission in the offered features remaining? How much stuff would that be? Does including a representative one of them cover most use cases? Would it in principle sound reasonable to include everything from the minimal set in the module?" This consideration seems to apply at a higher modular level. We propose only one or two new functions to be added. - "Is what you are adding something that is well understood in the literature? If not, how sure are you that it will turn out well? Does the method perform well compared to other similar ones?" Our article explains the basis of the method in elementary techniques. It does not require a sophisticated theoretical framework. Compare Theorem 6 from the article with the meat of the implementation near line 239 here: https://github.com/nadeemlab/coincidencetest/blob/main/coincidencetest/_coi… As for comparative performance, we have added support for the brute-force computational strategy by way of direct demonstration that the test is, indeed, an exact test in the standard sense. Our closed formula of course performs faster. If anyone here expresses interest in reviewing this proposal (the final criterion, "Someone wants to review it"), I will open a PR and get the ball rolling. As usual, please also feel free to voice any objections. -Jimmy

2 3