Add maxima interface to the step matrix finder (557d7ed9) · Commits · Andreas Klöckner / leap

@inducer Thanks! I have this capability in my branch and will try feeding it to Maxima

@inducer is this supposed to be returning the step matrix? Right now, what it's giving me is a vector of initial values and a matrix that isn't square - I'm not sure if the latter is supposed to be the step matrix, or if it's supposed to be something else altogether.

What I thought this would produce would be a step matrix I could let Maxima crunch on, but maybe there's something I've misunderstood.

Cory Mikida inform@tiker.net writes:

@inducer is this supposed to be returning the step matrix? Right now, what it's giving me is a vector of initial values and a matrix that isn't square - I'm not sure if the latter is supposed to be the step matrix, or if it's supposed to be something else altogether.

What I thought this would produce would be a step matrix I could let Maxima crunch on, but maybe there's something I've misunderstood.

This gives you the expression for one step. I.e. starting from initial_values, what would the state vector be after one step. The step matrix is the Jacobian of that vector with respect to initial_values. Taking those derivatives is exactly the step that leap's built-in symbolic stuff was having trouble with, and that's why that's getting outsourced to Maxima.

@inducer Thanks for the clarification on this. As an update - I'm now able to use this output in Maxima to generate a symbolic step matrix no sweat, and it seems to accomplish this more quickly than Leap would. However - unless I'm not using Maxima correctly - it appears to struggle greatly with finding the eigenvalues of that matrix at specific timesteps, as I'll need it to do in order to characterize the stability bounds...we can talk more about this tomorrow.

Heh, yeah. What Maxima tries to do is find eigenvalues symbolically, which is massively harder than trying to do so numerically (and frequently just outright impossible):

For example:

(%i1) A:matrix([sqrt(2), -1], [1, sqrt(2)]);
                             [ sqrt(2)    - 1   ]
(%o1)                        [                  ]
                             [    1     sqrt(2) ]
(%i2) eigenvalues(A);
(%o2)               [[sqrt(2) - %i, %i + sqrt(2)], [1, 1]]

That's a much cleaner result than "somewhere near 1.4142 -1.0001j, maybe, if the floating point gods were in a good mood". :)

This thread:

http://maxima-discuss.narkive.com/k5PoEp0k/compute-eigenvalues-numerically

seems to indicate that you can

load("lapack");

and then call zgeev(matrix). (Why zgeev? Because that's how LAPACK says double precision complex (z) general-format matrix (ge) give me the eigenvalues/vectors (ev) in the six characters that F77 allows for function names.)

Update: that worked for me. load("lapack") goes away for a while and does a bunch of scary-looking compilation stuff. But once it was done, stuff worked.

@inducer I have this working in Maxima. I'm able to compute the eigenvalues by loading Lapack routines, but it's still doing so very, very slowly (still taking too long for me to make use of this code). Unless you have any immediate thoughts on that, we can talk about it on Thursday, just thought I'd give you an update in the meantime.

What's the bit that takes long? zgeev? How big is the matrix?

(zgeev should take a couple seconds at most for anything that's smaller than 1000x1000.)

(I.e. same as Matlab eig--it's the same routine.)

The matrix I'm currently running with is 576 x 576 - right now I just have zgeev stuck in a while loop that tests an increasing timestep until an eigenvalue crosses 1, then it exits. By my estimation, each run through the loop (which is just a substitution into the step matrix followed by a zgeev call) is taking about 30 seconds.

I'm aware that my search algorithm isn't particularly sophisticated, so that's probably the next thing I'll try to tune up, but either way, if each eigenvalue solve is taking that long there's only so much I'll be able to do, especially since this is still for a step ratio of 2 and the matrix will continue to grow as that ratio increases.

Is it the substitution or the zgeev that's slow? (I'm guessing it's the substitution.)

The thing is that once we have the matrix, we could come full-cobbled-together-circle and reimport the matrix into pymbolic, which can then generate Python code to evaluate the matrix, which might be a good bit faster. Bonus: You could reuse the search algorithms that are already in leap.

It actually looks like the zgeev is the slow part - I tried taking the substitution part out of my loop and it's still taking the same (really long) amount of time.

I'd be all for reimporting the matrix in order to use the code I'm more familiar with - I'm just not sure what that would involve coding-wise on your end. I think at this point I'd be more comfortable showing you the very short maxima code I have to make sure I'm not doing anything boneheaded before trying to stick this thing back into the Python, but then again for all I know reimporting the matrix wouldn't involve much work at all.

That sounds like good news, because it should be fixable. Let's chat tomorrow.

@inducer An update on this/a couple of questions:

I'm trying to use the code you implemented above to get the Maxima expressions needed to get the step matrix, then feed them to Pymbolic's Maxima interop stuff, but when I try to parse the thing that Leap tries to give Maxima, it complains about it being too long - there appears to be a character limit in the parse of 2048, which we're far surpassing. I could try to figure out how to split the two assignments that the Leap code returns into separate strings and feed them in one by one, but I doubt that'll work because they'll still be too long. Is there a particular motivation for this line length limit? I vaguely recall you talking about some sort of memory limit past which this interop fails to work well - is that why this line limit length is implemented? Can I increase it?
I'm pretty sure I have a good grip on how to feed stuff to Maxima using Pymbolic's interop stuff, but it's not clear how to pass things in the other direction, i.e. from Maxima back to my Python code. Am I going to have to write my own code (i.e. add something to Leap) for this?

Cory Mikida inform@tiker.net writes:

@inducer An update on this/a couple of questions:

I'm trying to use the code you implemented above to get the Maxima expressions needed to get the step matrix, then feed them to Pymbolic's Maxima interop stuff, but when I try to parse the thing that Leap tries to give Maxima, it complains about it being too long - there appears to be a character limit in the parse of 2048, which we're far surpassing. I could try to figure out how to split the two assignments that the Leap code returns into separate strings and feed them in one by one, but I doubt that'll work because they'll still be too long. Is there a particular motivation for this line length limit? I vaguely recall you talking about some sort of memory limit past which this interop fails to work well - is that why this line limit length is implemented? Can I increase it?

Well, that line limit is currently hardcoded in pymbolic.interop.maxima, based on my observation of roughly where stuff started breaking. You'll have to remove that limit, submit bigger output, observe what breaks, and then fix that--in precisely that code.

I'm pretty sure I have a good grip on how to feed stuff to Maxima using Pymbolic's interop stuff, but it's not clear how to pass things in the other direction, i.e. from Maxima back to my Python code. Am I going to have to write my own code (i.e. add something to Leap) for this?

The reverse direction involves parsing Maxima output. An initial shot at this is in that same file, but I make no promises that it'll suit your purpose.

Happy to help if you run into trouble.