Optimizing time spent on numerically solving ODEs?

[ijustlovemath]

Thanks! Yeah, sorry about the missing classes; to get you a fully functional copy I'd need to share code I'm not allowed to give up. In response to the parsing, I will definitely try and eliminate duplicate efforts. One of the strange little caveats of this particular simulator is that it needs to be run just once at a time; instead of putting t = 0 and t = 10 with step h = 1, the end time has to be equal to start time + h, as without this you lose the "custom curve" information that's done in the F(X,t) VI. Things are being continually refit hence the RK4 simulator needs to be (unfortunately) continually called. Maybe I can figure out some inplace element tricks with that parsed data array to keep continual processing to a minimum.

WRT modifying the RK4 VI itself, what's the procedure for modifying things from standard libraries? Just make a copy and save it as a separate VI?

[johnsold]

Yes. When you start modifying vi.lib VIs, make a copy (preferably with a modified name) and put it somewhere else. Make it part of your project or a reuse library.

You never want to modify anything in vi.lib because the next upgrade or LV repair will restore everything to NI default and you will lose all the modifications. Also, some other developer (or you) who wants to use the original VI for a differetn project would not have it available.

Lynn

[ijustlovemath]

Okay, so I ended up taking a look at this and man, that RK4 code is a mess!

I did some profiling and with the default rk4 solver, the formula parsing takes about 10x the time of a single iteration timestep. Because I'm running this single-step, the continuous parsing was taking a lot of the time. Because the formula is not fixed in time, eg there are "constants" in the formula that are actually time variant, I thought maybe I could just do the main data structure building during initialization, and do some inplace replacement as data values change.

I don't know if you're familiar with the RK4 VIs, but I looked into "Runge Kutta Check.vi" and it seems like the "Table" output is constant if you format the inputs properly. I've attached an lvproj where I did this.

The one issue is that if you change the sign of any of the controls, except for i2 (which is only there to affect the value of "alpha" and not its sign), the "Table" structure changes; I'll get into this a bit later. I think the "Table" uses some sort of reverse polish notation that's built iteratively during the parsing stage.

There's also the "results" array, which encodes the values of all the "storage" constants, as seen in "Compiler Tree Structure.vi". This uses the tree structure built during the previous stage of parsing to extract everything it sees as "not a variable" and store it in this array at the appropriate index. It only actually stores the magnitude of this constant as negative numbers are built using the negation operation in the "Table" rather than just kept with the numbers in the "results" array.

It's not 100% clear to me right now why certain things go to certain places, as changing the sign of an element will swap around where things are. Luckily it makes no distinction between 0 and positive numbers, so a number becoming 0 is not important. Try changing "Ib" from 10 to -9 and you'll see the 10 in the second and third rows turn into a 9 and get shuffled around.

Now, all of my interpolation functions are actually going to have strictly non-negative results, so the constants are not going to change sign during execution.

So my approach is to note the indices where a certain constant is stored, and update them as they change during the simulation via an inplace element structure.

If you think you have some inkling as to why the elements get shuffled around and array sizes change after swapping signs, or know anyone who might, that would be very helpful in making the end code more robust.

Some alternate approaches:

1. When things do change sign, just reparse.

2. Add all of my time variant functions in as variables whose derivative is 0, eg instead of (g,x,i1,i2), have (g,x,i1,i2,i3,i4,v) and just set F(X,t) = 0 for the additional variables, as I will set them up via X0.

Thoughts?

[johnsold]

Now you can see why I did not try to put together a complete solution! Reverse engineering that thing is a nightmare.

I think you are on the right track. Even a modest amount of reduction of parsing may speed your process enough to be worthwhile.

I will look at your code tomorrow. 

Alternative 2 is something I considered but was not certain enough of your equations to be able to evaluate it. 

Lynn

[johnsold]

I think I found a way for your VI to extract the indices of the elements in the results array.

The Indices2 array is a 2D array of cluster of 1D array of cluster of points [m,n]. Each point [m,n] represents the indices of an element in results with the same value as element [i,j]. Indices2 is the same size as results. Every element of Indices2 contains at least one point because it always finds itself. The 79 zero elements in results result in the array at those elements having 79 points.  This is redundant but the array is not very large and this method is easier than a lot of fancy bookkeeping to reduce the redundancy.

Pseudocode:

For each i,j

   loop over all m,n

          if results[m,n] = results[i,j], 

               then append [m,n] to temporary indices array

  end m,n loop

  insert temporary indices array into Indices2[i,j]

Lynn

[ijustlovemath · Accepted solution]

Thanks, that was really bugging me! So I ended up dividing down about as deep as I was willing to go to optimize the RK4 VIs for "single step" simulation.

Here's what I found:

1. Overall results, one "cycle" of simulation ised to average 46ms, now it averages 5ms. Performance boost of 9200%, which is pretty awesome.

2. Here's how the main RK4 solver looks now:

3. Inside the "Iteration Setup" vi:

3. Inside the F(X,t) vi:

This is the most different from the original, which you posted previously. Per your advice, I do all the "code" conversion on initialization, and just use an accessor to pull the data from the class. Similarly, according to my other discoveries, I computed the "Table" and "Result" data structures on initialization and just stored them in the class private data for later access.

4. The "Phi" calculation:

I tested three different methods for this one: autoindexed formula nodes, autoindexed LV primitives, and Array-based LV primitives. The order from fastest to slowest was: Array-based primitives, autoindexed primitives, Formula node, each with a ~100x speed increase for large-ish inputs, and about a 10x for smaller inputs, so lots of time saved here.

5. Here's where we actually execute the parser:

There may yet be a little bit of duplication to reduce here, but it only takes a few ms and only runs once so I'm not too worried about it. I may also take all the controls outside of the case structure as I hear that has some performance increases. The data structure it computes gets stored in the class private data.

6. The "result" accessor (finally a use for dynamic accessors that I can justify in my data model!) for my ODE model specifically has this internally:

Note this gets called in each call of F(X,t), so another optimization may be to have an "update" that involves doing this calculation at the end of each simulation iteration. I also check what the model has stored, which isn't updated until the iteration is complete, and hence there may be some inconsistencies in the k_i of the RK4 approximation. I intend to more throughly review this today. A better approach may be to have this accessor take in the (t, X) array to allow my model and others like it to actually inspect what X is to make the given replacements in the "results" data structure. I'm also unsure if, in the future where I may have multiple indices in which to replace the constant in the data structure, I should be using an inplace Element structure inside some sort of autoindexed for loop, where the autoindexing thing is the 2D array of indices you helped me find earlier. That seems like it may be bad practice, though I'm not experienced anough to know for sure.

7. This actually really simplified my simulator code, which is an added plus. This is mostly due to the fact that I decided to use approach 2 of my previous post, which made for me having to do a hell of a lot fewer inplace replacements in the "results" data structure.

Anyways, thanks for the advice and the nudges in the right direction. I'm going to mark this as solved.

[johnsold]

Very nicely thought out and explained.

I saw in one of your other posts that you are relatively new to LV, so this is quite impressive.

I hope that someone at NI will bring this thread to the attention of the developers who work on the math VIs. Considering the iterative nature of ODE solvers, the improvments you have made should be considered for incorporation into vi.lib. The original VIs probably will be kept for compatibility with legacy code and for places where the parsing is not an issue.

Lynn

[ijustlovemath]

Thanks, that means a lot! I've been using LV almost daily for just over two years and recently got my CLAD, so it's nice to see that my studies are paying off!

I think there are different use cases here; the general case which currently exists is there because most people will plop a solver down, put in a start and end time and simulate for ~1k iterations, at which point the time the parser takes is pretty negligible. Most people also don't have a set of equations that changes in time either, eg max(0,I2), so the iterative nature of reparsing isn't there, although for people in my situation, this technique should help loads. When I was going deep into the RK4 code I noticed that their error checking involved unbundling the status element into a normal case structure, so this code probably is very very old. I think that in the future they could change the parser into a functional global that stores the parsed parameters and only reparses if F(X,t) or X have changed. You'd get the benefits of avoiding continuous reparsing without having to involve OOP. I think better documentation on the data structures involved would be a benefit to anyone looking to make their code more lean as well. There's also a few pretty obvious things that never should have made it past testing, like the auto indexed formula node, but I guess that's up to the maintainers of vi.lib. At the very least, it's open to the very bottom VIs, which is more than you can say for most (cough, Database Connectivity Toolkit).