Sorry to pick up on this old thread, but I too was looking for a faster expm implementation, and had some joy. Thought I'd share...
The expokit library seems to be a highly regarded matrix exponentiation library (that Matlab probably uses), which is written in Fortran. I downloaded it a while ago and compiled it with f2py; first time I've used f2py and it worked pretty much out of the box!
The Expokit source code is available from:- http://www.maths.uq.edu.au/expokit/download.html
Expokit just needs to be compiled and linked against LAPACK and BLAS libraries, so I've got two versions: one linked against NVIDIA's libcublas and another linked against Intel's libmkl (most people would link against blas and lapack libraries from ATLAS, I guess). The library can be specified manually, but it's probably easiest to just let numpy choose (like the command below does).
I've had it kicking around for a while, but just checked, and a compile command that just worked:-
f2py --fcompiler=intelem -c expokit.f -m expokit --link-blas_opt --link-lapack_opt
Unless you use Intel's 64bit Fortran compiler, you'll want to remove or change the --fcompiler option.
This creates the Python shared module: expokit.so This assumes a working numpy installation, and fortran compiler.
I've only needed to use Padé approximation myself, which is done with the subroutines [ZD]GPADM. These should give equivalent results to scipy.linalg.expm, but they take a lot more function arguments.
time_expm.py, attached, has an example wrapper function and runs some basic timings against scipy.linalg.expm.
The only change I made to the Fortran code was to add f2py style comments, so the output variables do return changed, and can be seen when inspecting with the numpy.info function. e.g.
>>> import numpy as np >>> import expokit >>> np.info(expokit.dgpadm) dgpadm - Function signature: dgpadm(ideg,t,h,wsp,ipiv,iexph,ns,iflag,[m,ldh,lwsp]) Required arguments: ideg : input int t : input float h : input rank-2 array('d') with bounds (ldh,m) wsp : input rank-1 array('d') with bounds (lwsp) ipiv : input rank-1 array('i') with bounds (m) iexph : in/output rank-0 array(int,'i') ns : in/output rank-0 array(int,'i') iflag : in/output rank-0 array(int,'i') Optional arguments: m := shape(h,1) input int ldh := shape(h,0) input int lwsp := len(wsp) in/output rank-0 array(int,'i')
For further information on the subroutine names and arguments, see the Expokit download link, above, which describes each subroutine, and provides separate links to each subroutine's source code. The Fortran source code is where to find the best documentation of each function and its arguments.
In terms of merging this into scipy, some work would need to be done...
1) Wrapper functions would need to be created to conform to scipy's expm API; time_expm.py, attached, demonstrates python wrappers for dgpadm and zgpadm. Using PyFort to automatially write and compile C wrappers should give further, marginal performance improvements. 2) Unit-tests on the wrapper functions. I think scipy's expm unit-tests could be used here. 3) The scipy build systems would need to be informed about it. 4) scipy.linalg.__init__.py should probably try and import expokit wrapper functions over the expm, expm2 and expm3 functions. If that fails, fall back to the original scipy versions.
Performance improvements? These are the results I got when just running it, averaging the times over 3 runs. Note, both scipy and expokit use the same BLAS dot product functions, which does a lot of the heavy lifting; over multiple processor cores, I might add.
----- Mean time +- S.D. ----- Array size Expokit scipy..expm 25 0.149 += 0.008 0.410 += 0.025 50 0.341 += 0.011 0.678 += 0.019 75 0.650 += 0.049 1.282 += 0.030 100 1.372 += 0.087 1.995 += 0.032 125 1.972 += 0.029 3.007 += 0.145 150 6.038 += 0.062 7.933 += 0.217 175 9.169 += 0.050 11.785 += 0.593 200 12.049 += 0.190 16.281 += 0.293
Hope this might be of help to someone. Kind regards, Alex _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
Alex this is some very nice work. Unfortunately I don't think the
Expokit licensed (http://www.maths.uq.edu.au/expokit/copyright) is
compatible with the BSD license SciPy uses. Specifically approval
must be obtain before commercial use. Including the package in the
Topical Software list and/or the SciPy Cookbook/SciPy Central might
be prudent.
- Jonathan Helmus On 11/07/2012 09:47 AM, Alex Leach wrote:
_______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
On Wednesday 07 Nov 2012 10:04:55 Jonathan Helmus wrote: > Alex this is some very nice work. Unfortunately I don't think the > Expokit licensed (http://www.maths.uq.edu.au/expokit/copyright) is > compatible with the BSD license SciPy uses. Specifically approval must > be obtain before commercial use. Including the package in the Topical > Software list and/or the SciPy Cookbook/SciPy Central might be prudent. > > - Jonathan Helmus
(sorry, replied directly before)
Hi Jonathan,
Good spot about the license issue. That's a shame, but seems like the author might give SciPy permission to use it, if asked..
Including it as a package, documented build method, or otherwise would be cool. If someone wants to go ahead and do that, that's fine; I don't think I have the access rights (or the time, really) to do any of that myself.
One other thing; just did a diff to see if I'd made any other changes. Looks like I also had to edit some calls to the matvec functions. Converted C code fails during compilation without them. Universal diff attached.
Cheers, Alex _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user expokit.udiff (11K) Download Attachment |
On 11/7/2012 10:57 AM, Alex Leach wrote:
> seems like the author might give SciPy permission to use it, if asked. So, will you ask? Just fyi, I have had good luck with such requests, when I ask if BSD/MIT would be possible and explain why it is needed for SciPy. Alan Isaac _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Bugzilla from beamesleach@gmail.com
On 7 Nov 2012 15:58, "Alex Leach" <[hidden email]> wrote: Note that just asking for SciPy to have permission is insufficient. We actually need it to allowed for use under bsd-like terms. This is because we need to make sure that everyone who uses scipy also can use it. The current license makes it impossible to use expokit and GPL code in the same program, eg. -n > _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
On Thursday 08 Nov 2012 07:19:51 Nathaniel Smith wrote: > > Note that just asking for SciPy to have permission is insufficient. We > actually need it to allowed for use under bsd-like terms. This is because > we need to make sure that everyone who uses scipy also can use it. The > current license makes it impossible to use expokit and GPL code in the same > program, eg. >
Kind email sent. Fingers crossed.
Alex _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Nathaniel Smith
Great news guys. The Expokit author, Professor Roger Sidje, is perfectly happy for Expokit to be used in SciPy. He just likes to keep a tab on things, for appraisals and such like.
His full response, and my original message, is below:-
Alex,
Thanks for your interest.
Please note that the wording of the license is merely a way to keep track of where it is used rather than an attempt to restrict its usage.
The license is simple -- any other open/closed source can use the software any way they want.
All the license does is to allow me to know/have a bullet point in those end-of-year activity reports for the powers that be, telling that the software is found useful.
I have had similar requests in the past re:GPL. All I say is, feel free to use/modify the software any way you want. This e-mail of yours is exactly the kind of expression-of-interest indicated (even for use in a commercial setting).
Because I have moved from UQ (in Australia) to UA (in the USA), making changes to Expokit won't be convenient for a while. Please consider this response as a plain acknowledgment that Expokit can be embedded/released under GPL. You can add the following in the meantime:
"Alternatively, EXPOKIT may be used under the terms of the GPL."
By so doing, it would mean that you have chosen to branch on the GPL-alternative, while allowing others to do what they want on the other alternatives (and it is understood that, per GPL, if someone ever gets the GPL alternative or its derivatives, they will abide by GPL).
In particular, possible future versions of Expokit not derived from your copy will remain free to be licensed any way (which is an option that commercial packages want -- even though they use EXPOKIT without royalties).
Thanks --- RBS
On Thu, Nov 8, 2012 at 5:44 AM, Alex Leach <[hidden email]> wrote: > Dear Professor Sidje, > > I hope this email finds you well. > > As way of introduction, I am a PhD student studying Computational Biology at > the University of York. I have recently used Expokit for modelling > polymerisation in a specific chemical reaction, and had some good > experiences with it. > > I thought I'd share my joy with the SciPy community, who have implemented > similar algorithms for computing matrix exponentials, but in Python. After > writing a small wrapper function (and making some very minor changes to the > Expokit source code), Expokit's dgpadm worked as a drop-in replacement for > their `expm` implementation, but was about twice as fast. > > The exact changes I made can be seen in expokit.udiff, which I sent as an > attachment with this email to the SciPy list. (The SciPy mailing system made > the nasty file name: attachment-0001.bin.) > > Impressed with the speed improvements, the SciPy community would be > interested in using Expokit in precedence of their original implementation. > > However, there is an issue with licensing... Python is GPL source code, and > SciPy uses a BSD license. I'm sure you know more about these licenses than I > do, but from what I gather, the Expokit Copyright is restrictive for > corporations who would wish to use the software for commercial purposes, > whereas the BSD license is non-restrictive in this sense. > > So, as you've probably guessed, I'm writing to ask your kind permission, on > behalf of the SciPy community, to include and distribute the Expokit source > code, as a part of the SciPy source code distribution. > > With your permission, expokit.f would be included in future SciPy releases > and built as a Python shared module, which is compiled from intermediate C > code that includes the relevant Python headers, etc. The first link in this > message, above, shows the performance improvements I achieved, build > instructions and python test module, detailing how I did this. > > If you have any questions or requirements with regard to the usage of > expokit, please don't hesitate to ask. > > There would be some work to merge Expokit into SciPy, but there seems to be > a lot of interest in doing this, so if you give the green light, the > community would get this done in no time, I reckon. It's impossible to tell > how many tools would benefit from this change downstream, but I know that > some popular protein visualisation tools use SciPy internally, so I think > this would benefit a massive user-base. > > I look forward to hearing back from you. > Kind regards, > Alex > > > -- > Alex Leach BSc. MRes. > Department of Biology > University of York > York YO10 5DD > United Kingdom > www: http://bioltfws1.york.ac.uk/~albl500 > EMAIL DISCLAIMER: http://www.york.ac.uk/docs/disclaimer/email.htm _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
On 8 November 2012 23:19, Alex Leach - Gmail <[hidden email]> wrote: Great news guys. The Expokit author, Professor Roger Sidje, is perfectly happy for Expokit to be used in SciPy. He just likes to keep a tab on things, for appraisals and such like. It certainly sounds like he's happy to see the code reused, but I'm not convinced that his e-mail technically provides the necessary permission: > All I say is, feel free to use/modify the software any way you want. This statement appears to grant a very permissive license, but it doesn't mention any right to redistribute the code, which is key. > Please consider this response as a plain acknowledgment that Expokit can be embedded/released under GPL. Unfortunately, we can't integrate GPL licensed code into Scipy, as Scipy is BSD licensed. To integrate that code, we'd need to get a similarly permissive license. I've just found this handy page: http://www.scipy.org/License_Compatibility He might not be willing to give us that - he wants to ensure he knows about who's using his software, which is a requirement you can't embed in a permissive license. I think it's worth a try, though. A couple of points you could use: - We're also keen to know who's using our software, although we don't make it compulsory for commercial users to tell us. E.g. IPython maintains a list of projects using IPython, which we use to demonstrate impact. - He could grant a BSD license to just the bit of the code we're interested in, rather than all of Expokit. Kudos for working on this, Thomas _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
On 9 Nov 2012, at 00:57, Thomas Kluyver <[hidden email]> wrote:
All good. Just got confirmation back; below. On 9th Nov 2012, at 14:44, Roger B. Sidje wrote: As I indicated the license is very permissive indeed. Just substitute GPL with BSD in my earlier post. You can branch with your BSD copy, make change and redistribute SciPy with those changes any way you guys want, while allowing others to do what they want with the other alternatives not derived from your branch. --- RBS _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
On 9 November 2012 15:12, Alex Leach <[hidden email]> wrote:
That sounds OK to me :-) Thomas _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
Free forum by Nabble | Edit this page |