Numerical solver design is a topic that I know next to nothing about. I always say I know a tiny bit of MATLAB very well, but lately in my push to become MATLAB-independent, I’ve been trying to remove bits of proprietary code from my own simulations. (My transition from MATLAB to Octave failed in part because of memory issues in Octave. So I’m now looking at Python as my potential free software savior.)
The biggest hurdle has been my usage of the suite of numerical solvers for ordinary differential equations (ODEs) in MATLAB. These are fast, vectorizable to some extent, and take advantage of some parallel processing on multicore machines. However, using their solvers presents a number of problems when one wants to understand the ins and outs of their algorithms. It’s not that this information isn’t there (using ‘type’ or ‘help’), but the problem with going through code that isn’t yours is an ugly one that isn’t going away soon. Additionally, it’s never been completely clear under what circumstances MATLAB’s solvers start to make use of multiple cores and why that process should be accurate, despite seeing evidence of multicore usage at some level of scale in my own code, with ode23s. And there are some things that I don’t know if MATLAB’s solvers can do well. For instance, I have MxNxT matrices of a dynamical variable connecting M model cells of one type to N model cells of another type, saved over all T time steps. Tracking all of these variables and the some ~20 other variables of varying lengths becomes a real code nightmare in vector forms without name labels, which are available in structs. Of course, there is speed in vectors that isn’t available in structs — speed that I am willing to sacrifice for sanity.
The real bottom line for me is being able to write some MATLAB code that I can easily translate into a Python environment eventually. This process will take a lot longer than I had hoped. By writing solvers now, I hope to be able to take them with me and get some insight into solver design for future projects.
Using structs in a solver. In MATLAB’s ODE solvers, they require a strict type of vectorization that has time in one dimension and all of the dynamical variables in the other dimension. I have synaptic weight matrices that are MxNxT, for M of cell A going to N of cell B, throughout time T. Because of the dimensions of these are tricky to code in MATLAB directly without a lot of transposes and other confusing situations, I was hoping to use structs to name my matrices and make it easy to write vector field equations in a clear manner without having to write confusing converters. It’s not completely out of laziness that I want to avoid doing this; rather, another conversion step introduces one more place where error can be introduced.
It turns out that using structs is not a good idea with the separate solver/vector field design, since it seems like you’d spend too much computational time translating anyway. Translating once in a difficult way but gaining speed is far better than translating thousands of times in an “easy” way, at some unknown expense of speed. I might try next a two-matrix approach: one for class A cells and one for class B cells, which does not make this a general solver but will allow for easier vectorization.
Variable passing to a function using a solver. In Matlab, one thing I had trouble learning early on was the method by which one passed variables to the vector field file via the solver. The syntax might look like:
[t, y] = solver('vf', t_range, y0, , var1, var2, ...);
For the outputs, we have simply the t vector and the numerical solution y for the system, which is specified in the vector field function, in this case ‘vf’. The range of times and the initial condition y0 should be passed to the solver, which uses ‘vf’ to solve the system. Now if the function ‘vf’ requires inputs, which is not uncommon, in order to pass them to ‘vf’, you have to pass them through the solver first. If you want to pass ‘var1’ and ‘var2’ to ‘vf’, notice that you cannot simply specify these variables after the ‘y0’. Instead, in this example, there are empty brackets. These empty brackets are actually part of the solver’s input syntax to allow options to be passed to the SOLVER and not to the VF. The reasons for this are very clear to me, as I code my own solver and recognize that this is exactly how my solver ended up being coded. The solver needs the ability to accept options, and they need to come into the fixed location here. Any additional arguments passed through the solver can then be understood as variables to the vector field. While it seems reasonably obvious, I understand it much better after having tried to figure out this design from first principles.