Mailman 3 April 2016 - SciPy-Dev

scikit-kinematics
by Thomas Haslwanter April 9, 2016

April 9, 2016

Hi, I am thinking about taking all the kinematics-related modules from my library thLib (http://work.thaslwanter.at/thLib/html/#modules), and release it as a new scikit: "scikit-kinematics". Currently there are modules on: - quaternions - rotation matrices - vector calculations - IMU (inertial measurement units) - markers (optical tracking systems) - viewer (for interactively investigating time series) Any opinion?

1 0

thLib - 3D Kinematics and Signal Processing
by Thomas Haslwanter April 7, 2016

April 7, 2016

Hi, over the last few years I have developed a Python library, primarily aimed at 3D kinematics, signal processing, and miscellaneous small items that I use regularly, like a waitbar. It is currently available at https://pypi.python.org/pypi?%3Aaction=pkg_edit&name=thLib Somewhat to my surprise, it has become pretty popular (by my personal standards), with almost 9'000 downloads last month. Since most parts of it are by now quite stable, I am thinking about how best to continue. I can 1. Make it a github-project, and see who is interested in contributing. 2. Improve the testing, and release a 1.0 version on my own. (My main limitation is that I don't have a Mac, and testing compatability with Macs is difficult for me.) 3. Turn it into a "scikit-kinematics" (or some other suitable "scikit"-name) and/or contribute the code to scipy, if others would find it suitable. I would appreciate feedback and opinions. Thanks, Thomas

2 2

Edge case bug in stats.triang?
by Andrew Nelson April 6, 2016

April 6, 2016

I've been investigating the stats.triang.pdf over it's domain with different modes: >>> from scipy.stats import triang >>> triang.pdf(0, 0) 2.0 >>> triang.pdf(1, 1) nan Shouldn't the last value be 2? -- _____________________________________ Dr. Andrew Nelson _____________________________________

2 1

Introducing scipy.optimization.tgo; a global optimization routine.
by Stefan Endres April 6, 2016

April 6, 2016

Hello all, I’d like to present this implementation of the Topographical Global Optimisation (TGO) method, a powerfull global optimisation method first proposed by Törn (1990) [1]. The implementation is written entirely in python using mostly boolean numpy arrays to implement the core algorithm and using scipy.optimize.minimize to find local minima. https://github.com/Stefan-Endres/tgo TGO is a clustering method that uses graph theory to generate good starting points for local search methods from points distributed uniformly in the interior of the feasible set. These points are currently generated using the Sobol (1967) [3] sequence. The semi-empirical correlation by Hendorson et. al. (2015) [2] for k integer defining the k-t matrix. The latter publication used an (unpublished as far as I can tell) C/C++ implementation which was comparatively succesful against commercial software (MIDACO 4.0). More details can be found in the docstrings in https://github.com/Stefan-Endres/tgo Some features in brief. - Global optimisation of multivariate functions with - Bounded or unbounded variables. - Optional non-linear constraints. - Returns an ordered list of all local minima with the res.xl attribute and corresponding function values res.fun (this is something that I’ve found very useful in my own research). - Highly customisable to increase performance on difficult optimisation problems. - The functions were written entirely in python based on algorithms available in open literature [2] As a quick visual illustration of the algorithm’s capabilities, a few results for common optimisation problems and their results: http://imgur.com/a/'ll MS0hG <http://imgur.com/a/MS0hG> I’ve also run a comparison of tgo with the other two global optimisation currently available. http://pastebin.com/8K72Xve8 Note in all cases tgo found the correct answer, details on my system can be found at the end of the pastebin (it’s a terminal session dump). Also I apologise I’m not very familiar with the basinhopping algorithm so I didn’t know how many iterations where optimal for the problems and stuck with the default. Details on the functions used and their references can be found in the docstrings in the unittest script: https://github.com/Stefan-Endres/tgo/blob/master/test__tgo.py Potential issues The most important barrier I forsee is the large file size of the direction numbers used in generating the Sobol sampling points (new-joe-kuo-6.21201 <https://github.com/Stefan-Endres/tgo/blob/master/new-joe-kuo-6.21201>) at 1.8MiB. The code to generate the Sobol points itself is not an issue, it’s a small script that was translated to Python by one of my advisers for the project (repository: https://bitbucket.org/alchemyst/sobol/src) from the original C++ code under a scipy compatible BSD style liscence ( http://web.maths.unsw.edu.au/~fkuo/sobol/). I suggest the following solutions to this: - Use latin hypercube sampling (LHS) to generate the initial points instead. LHS is already used in _differentialevolution.py (I’ve tried to do this, but I’m not sure how without initiating the entire DifferentialEvolutionSolver object). - Trim the data file to only including Sobol sequence sampling for problems of up to 1000 dimensions. It’s easy to write something to download higher dimensional data when a user specifies problems higher than 1000 dimensions. I think the sobol_points function is more appropriate numpy.random in the first from where it can be imported to tgo, but I'm not sure if it's wanted there. Everything else was written from scratch so there shouldn’t be any licensing issues. I’ve tried to follow the guidelines given in http://docs.scipy.org/doc/scipy/reference/hacking.html as close as possible with this repository: https://github.com/Stefan-Endres/scipy/branches It’s almost ready for a PR, in particular I don’t really know where/how to do the following checklists. - Is the new functionality tagged with .. versionadded:: X.Y.Z (with X.Y.Z the version number of the next release - can be found in setup.py)? - Is the new functionality mentioned in the release notes of the next release? - Is the new functionality added to the reference guide? References 1. Törn, A (1990) “Topographical global optimization”, Reports on Computer Science and Mathematics Ser. A, No 199, 8p. Abo Akademi University, Sweden) 2. Henderson, N et. al (2015) “A new look at the topographical global optimization method and its application to the phase stability analysis of mixtures”, Chemical Engineering Science, 127, 151-174 3. Sobol, IM (1967) “The distribution of points in a cube and the approximate evaluation of integrals. USSR Comput. Math. Math. Phys. 7, 86-112 If you wish to learn more about the method I highly recommend Hendorson et. al.’s (2015) article which presents the material in an introductory fashion. Please let me know what you guys think! I will keep working on this if you think it’s something that could be included in scipy. Thank you for reading. Regards, - Stefan Endres Postgraduate Student: Institute of Applied Materials Department of Chemical Engineering, University of Pretoria Cell: +27 (0) 82 972 42 89 E-mail: Stefan.C.Endres(a)gmail.com St. Number: 11004968

6 11

scipy.sparse.linalg.spsolve surprising behaviour for large sparse matrices on Linux systems
by Alain Batailly April 4, 2016

April 4, 2016

Hello, I am computing the solution of a linear system Ax=b with A a large (typically 200,000 lines and columns for the associated dense matrix) sparse matrix and b a sparse matrix of about 100 columns. When I run my code on Windows systems (Python *2.7*, scipy *0.14.0*), the following command from scipy.sparse.linalg import spsolve...Temp = spsolve(A.tocsc(),b.tocsc()) runs smoothly and requires about 7 GB of ram memory. Running the *exact same code* with the *exact same matrices* on Linux systems (same CPU, same amount of RAM memory: 64 GB, Linux Mint *17.3*, python *2.7*, scipy *0.13.3*) requires more than 20 GB of ram memory and it crashes with the following error message: <function umfpack_di_numeric at ...> failed with UMFPACK_ERROR_out_of_memory (see 1 <http://i.stack.imgur.com/ezdPi.png>) Because this error is os dependent, I ruled out any issue regarding the matrices *A* and *b*, and I am trying to find a fix specific to Linux... But I do not know where to start... Does anyone would have any idea of what is going on ? And why would such a problem be specific to Linux systems ? After further investigations, it seems that the latest version of BLAS on Linux Mint is quite old: 1.8.2. On Windows, I use BLAS 1.9.1. When using the test_numpy.py file available here: https://gist.github.com/osdf/3842524#file-test_numpy-py I notice very significant differences between Linux and Windows: *Linux* : version 1.8.2, maxint 9223372036854775807, dot: 0.76 s - *Windows* : version 1.9.1, maxint 2147483647, dot: 0,037 s. I realized that the problem may be hardware related. Indeed, an older PC, with the exact same libraries on the same Linux Mint distribution (Rosa 17.3) provides much more satisfying results. The benchmark mentioned in the first update gives on this old PC: *Linux* : version 1.8.2, maxint 9223372036854775807, dot: 0,054 s. Would anybody have observed such behaviour ? Thank you. Alain

1 0

Re: [SciPy-Dev] SciPy-Dev Digest, Vol 150, Issue 4
by Stefan Endres April 4, 2016

April 4, 2016

>1. sample points from all possible solutions (e.g. sobol points or a simple grid) >2. discard points that don't look promising >3. perform local optimization on the remaining points. >And I guess the graph theoretic / clustering cleverness comes in step 2. > (Something along the lines of "keep only points where the K nearest neighbors all have larger function values", if i'm not mistaken.) Exactly right, in brief the topograph is constructed from an ordered list based on the euclidean distance between the points, the algorithm compares the function values of each point and assigns either a positive or negative sign depending on whether or not the it's higher or lower than the first point (implemented here as True booleans for positive signs and False boolean for negative flags). "Minimiser" points are now defined as any row in the topograph with all positive entries (ie all other points relatively "near" that point is above it). Of course, because the sampling is (pseudo-)random you don't want to compare all the points and have only one minimser, since that might not be near the actual global minima unless the objective function is convex. That's where the importance of the k_t integer comes in since it essentially determines the number points you compare and the number of minimisers you find (which are ideally not near each other). I've attached a brief presentation which outlines this construction more formally for those who are interested (this is cut from a phase equilibrium lit. review of please ignore the irrelevant references the material balance constraints etc.). >probably doesn't scale to very high dimensions. I don't believe so no, using a very high amount of sampling points would involve N x N distance calculations organised into an array of the same dimensions. The highest I've tested is unittest test_t910 which is only 7 dimensions. The highest dimensional problem tested in Henderson's article is only 20 dimensions. >It's nice that the bounds implementation seems to be fairly robust. It does seem to work very well with `SLSQP` so far, I've been using it in my work for a while and have never had any issues, lets hope it holds for more complex problems. In the future I'd also like to add equality constraints to the tgo routine itself (currently it's only used in the local minimisation step) which would adapt the sampling in a similar way to the inequality adaption (see the also the figures in the attached presentation as a visual illustration; I believe there is room for improvement here by somehow shifting the points uniformly into the feasible set rather than just discarding them)). >People would probably find it useful to have a standalone implementation of sobal sequences. Maybe it would fit in scipy.spatial? I agree, it should definitely not be embedded in tgo in any case, my adviser said he is still planning to do a few code optimisations (this repository: https://bitbucket.org/alchemyst/sobol/src), I will speak to him about it once he's finished. >In general, I think the approach is nicely complementary with basinhopping and differential evolution. I think it could be a good fit for scipy. Excellent! I will start cleaning the branch in the meantime then untill we figure out where to put the Sobol function. On 4 April 2016 at 19:21, <scipy-dev-request(a)scipy.org> wrote: > Send SciPy-Dev mailing list submissions to > scipy-dev(a)scipy.org > > To subscribe or unsubscribe via the World Wide Web, visit > https://mail.scipy.org/mailman/listinfo/scipy-dev > or, via email, send a message with subject or body 'help' to > scipy-dev-request(a)scipy.org > > You can reach the person managing the list at > scipy-dev-owner(a)scipy.org > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of SciPy-Dev digest..." > > > Today's Topics: > > 1. Re: Introducing scipy.optimization.tgo; a global optimization > routine. (Jacob Stevenson) > > > ---------------------------------------------------------------------- > > Message: 1 > Date: Mon, 04 Apr 2016 17:21:09 +0000 > From: Jacob Stevenson <jstevenson131(a)gmail.com> > To: SciPy Developers List <scipy-dev(a)scipy.org> > Subject: Re: [SciPy-Dev] Introducing scipy.optimization.tgo; a global > optimization routine. > Message-ID: > <CAEEL3uxfJMJWtU7UD-K80BQO_cS-NV0R0FxVDoHXbq30= > 2mtgQ(a)mail.gmail.com> > Content-Type: text/plain; charset="utf-8" > > I agree that the algorithm looks promising (though I also don't have access > to the papers). If I understand correctly, a high level summary is > > 1. sample points from all possible solutions (e.g. sobol points or a simple > grid) > 2. discard points that don't look promising > 3. perform local optimization on the remaining points. > > And I guess the graph theoretic / clustering cleverness comes in step 2. > (Something along the lines of "keep only points where the K nearest > neighbors all have larger function values", if i'm not mistaken.) > > Here are a few thoughts I have, please chime in if I'm making bad > assumptions. > > It starts with a grid-like search, and as such probably doesn't scale to > very high dimensions. But I think that's fine, it would still be > applicable to many (most?) applications. > > It's nice that the bounds implementation seems to be fairly robust. That > is currently a weakness of basinhopping. > > I'm curious how the automatic selection of hyper parameters works. I guess > we'll find out in Andrea's benchmarks. > > People would probably find it useful to have a standalone implementation of > sobal sequences. Maybe it would fit in scipy.spatial? > > In general, I think the approach is nicely complementary with basinhopping > and differential evolution. I think it could be a good fit for scipy. > > Thanks for the contribution! > Jacob > > > On Mon, 4 Apr 2016 at 16:45 Andrea Gavana <andrea.gavana(a)gmail.com> wrote: > > > Hi Stefan, > > > > The algorithm seems very interesting. I haven't gone through the > > papers (mostly because I don't have access to those journals ?), so I > > can't really say I understand the theoretical basis of it... > > > > However, I believe it would be nice to put TGO through the grinder as I > > did in the past with many other algorithms: > > > > http://infinity77.net/global_optimization/multidimensional.html > > > > And take a look to the overall performances. I'm always on the lookout > for > > GO algorithms available in Python (and actually running on Windows as > > well...), so I'll try to find some time soon enough to put TGO as well on > > the tables and pictures. > > > > Thank you for sharing! > > > > Andrea. > > > > > > On Monday, 4 April 2016, Stefan Endres <stefan.c.endres(a)gmail.com> > wrote: > > > >> Hello all, > >> > >> I?d like to present this implementation of the Topographical Global > >> Optimisation (TGO) method, a powerfull global optimisation method first > >> proposed by T?rn (1990) [1]. The implementation is written entirely in > >> python using mostly boolean numpy arrays to implement the core algorithm > >> and using scipy.optimize.minimize to find local minima. > >> > >> https://github.com/Stefan-Endres/tgo > >> > >> TGO is a clustering method that uses graph theory to generate good > >> starting points for local search methods from points distributed > uniformly > >> in the interior of the feasible set. These points are currently > generated > >> using the Sobol (1967) [3] sequence. The semi-empirical correlation by > >> Hendorson et. al. (2015) [2] for k integer defining the k-t matrix. The > >> latter publication used an (unpublished as far as I can tell) C/C++ > >> implementation which was comparatively succesful against commercial > >> software (MIDACO 4.0). > >> > >> More details can be found in the docstrings in > >> https://github.com/Stefan-Endres/tgo > >> > >> Some features in brief. > >> > >> - Global optimisation of multivariate functions with > >> - Bounded or unbounded variables. > >> - Optional non-linear constraints. > >> - Returns an ordered list of all local minima with the res.xl > >> attribute and corresponding function values res.fun (this is > >> something that I?ve found very useful in my own research). > >> - Highly customisable to increase performance on difficult > >> optimisation problems. > >> - The functions were written entirely in python based on algorithms > >> available in open literature [2] > >> > >> As a quick visual illustration of the algorithm?s capabilities, a few > >> results for common optimisation problems and their results: > >> http://imgur.com/a/'ll MS0hG <http://imgur.com/a/MS0hG> > >> > >> I?ve also run a comparison of tgo with the other two global optimisation > >> currently available. > >> > >> http://pastebin.com/8K72Xve8 > >> > >> Note in all cases tgo found the correct answer, details on my system can > >> be found at the end of the pastebin (it?s a terminal session dump). > Also I > >> apologise I?m not very familiar with the basinhopping algorithm so I > didn?t > >> know how many iterations where optimal for the problems and stuck with > the > >> default. > >> > >> Details on the functions used and their references can be found in the > >> docstrings in the unittest script: > >> https://github.com/Stefan-Endres/tgo/blob/master/test__tgo.py > >> Potential issues > >> > >> The most important barrier I forsee is the large file size of the > >> direction numbers used in generating the Sobol sampling points ( > >> new-joe-kuo-6.21201 > >> <https://github.com/Stefan-Endres/tgo/blob/master/new-joe-kuo-6.21201>) > >> at 1.8MiB. The code to generate the Sobol points itself is not an issue, > >> it?s a small script that was translated to Python by one of my advisers > for > >> the project (repository: https://bitbucket.org/alchemyst/sobol/src) > from > >> the original C++ code under a scipy compatible BSD style liscence ( > >> http://web.maths.unsw.edu.au/~fkuo/sobol/). > >> > >> I suggest the following solutions to this: > >> > >> - Use latin hypercube sampling (LHS) to generate the initial points > >> instead. LHS is already used in _differentialevolution.py (I?ve tried > >> to do this, but I?m not sure how without initiating the entire > >> DifferentialEvolutionSolver object). > >> - Trim the data file to only including Sobol sequence sampling for > >> problems of up to 1000 dimensions. It?s easy to write something to > download > >> higher dimensional data when a user specifies problems higher than > 1000 > >> dimensions. > >> > >> I think the sobol_points function is more appropriate numpy.random in > >> the first from where it can be imported to tgo, but I'm not sure if it's > >> wanted there. > >> > >> Everything else was written from scratch so there shouldn?t be any > >> licensing issues. > >> > >> I?ve tried to follow the guidelines given in > >> http://docs.scipy.org/doc/scipy/reference/hacking.html as close as > >> possible with this repository: > >> https://github.com/Stefan-Endres/scipy/branches > >> > >> It?s almost ready for a PR, in particular I don?t really know where/how > >> to do the following checklists. > >> > >> - Is the new functionality tagged with .. versionadded:: X.Y.Z (with > >> X.Y.Z the version number of the next release - can be found in > setup.py)? > >> - Is the new functionality mentioned in the release notes of the next > >> release? > >> - Is the new functionality added to the reference guide? > >> > >> References > >> > >> 1. T?rn, A (1990) ?Topographical global optimization?, Reports on > >> Computer Science and Mathematics Ser. A, No 199, 8p. Abo Akademi > >> University, Sweden) > >> 2. Henderson, N et. al (2015) ?A new look at the topographical global > >> optimization method and its application to the phase stability > analysis of > >> mixtures?, Chemical Engineering Science, 127, 151-174 > >> 3. Sobol, IM (1967) ?The distribution of points in a cube and the > >> approximate evaluation of integrals. USSR Comput. Math. Math. Phys. > 7, > >> 86-112 > >> > >> > >> > >> > >> If you wish to learn more about the method I highly recommend Hendorson > >> et. al.?s (2015) article which presents the material in an introductory > >> fashion. > >> > >> Please let me know what you guys think! I will keep working on this if > >> you think it?s something that could be included in scipy. > >> > >> Thank you for reading. > >> > >> Regards, > >> - > >> Stefan Endres > >> Postgraduate Student: Institute of Applied Materials > >> Department of Chemical Engineering, University of Pretoria > >> Cell: +27 (0) 82 972 42 89 > >> E-mail: Stefan.C.Endres(a)gmail.com > >> St. Number: 11004968 > >> ? > >> > > > > > > -- > > > > _______________________________________________ > > SciPy-Dev mailing list > > SciPy-Dev(a)scipy.org > > https://mail.scipy.org/mailman/listinfo/scipy-dev > > >

1 0

Sampling from arbitrary distributions
by Andrew Nelson April 3, 2016

April 3, 2016

Hi all, I have a need to take samples from an arbitrary continuous statistical distribution. By arbitrary I mean it doesn't have an exact functional form. Rather, the distribution was created by sampling a physical process and then histogramming. I then wish to randomly sample from this distribution for a monte carlo simulation. Whilst I could (in this circumstance) sample with replacement from my original measurement I thought it might be quite cool to subclass scipy.stats.rv_continuous to do it, see the gist for further info. https://gist.github.com/andyfaff/9a7930f3d93a99f478cdadae4d1bb3b2 This way one gets lots of other benefits for free (e.g. the CDF, PPF, calculation of statistical moments, etc). I was wondering if this might be a useful thing to add to scipy.stats. I'm not overly familiar with the stats module so I thought I'd ask this again on the dev list (I asked this question on the user list a couple of days ago). The class does store state because it needs to do interpolation, etc. It was suggested that the "Histogram distribution would be the continuous analog to the arbitrary finite number of points discrete distribution, so might fit as well." Don't get too caught up by it working on histogrammed data, it should be straightforward to modify it to act for a distribution that has no analytic functional form, but where PDF(x) is known for many x. A. -- _____________________________________ Dr. Andrew Nelson _____________________________________

2 2