persistent variable precludes simbiology acceleration, how can I keep an array available in memory

In a simbiology model, I use a spline function in a matlab .m function to set the value of a variable parameter. For example, the assignment rule might be: glucose = forced_glucose(time) The function is called a lot, and I'd like to minimize overhead. The data for the spline is in an xlsx file, and I'd like the function to only read the file once, and then retain the data in memory until I clear the function. Using the following works very well when one is running the model and developing stuff and testing constructs:

function glu = forced_glu(time)
  persistent pp
  if isempty(pp)
  raw = xlsread('my data file.xlsx')
  pp = spline(raw(:,1),raw(:,2));
  end
  glu = ppval(pp,time);
  end

Unfortunately, when one uses sbiofit, a warning is flagged noting that the model could not be accelerated. Apparently this is a limitation in the accelerate function, and is not surmountable. The limitation IS noted in the documentation set for "accelerate" where one is warned not to use persistent variables in functions used by simbiology, but there is no warning in the docs for "persistent".

I've been goofing around a bit trying to see if the global construct would work. There were a few glitches in my approach. I thought that someone might have run into this. Is there a way to ensure that a function used in simbiology rate expressions has a way to initialize on the first call, and to keep a structure (the piecewise polynomial structure 'pp') in accessible memory while running the model? I really don't want to process date into splines every call (several calls per time point).

In MATLAB, I could make the function use the pp as both an input and an output argument. But would that work with SB repeat assignment rules? Can I say [simbio_parameter_value pp] = forced_glucose(time,pp), and have the output pp stored in matlab memory?

Ideas? Thanks!

Answer 1

[Arthur Goldsipe]

This is a timely question. The note about persistent variables being incompatible with acceleration is incorrect. (I think that used to be a limitation?) I do still encourage you to use them with caution, but pre-loading data for interpolation is a perfectly valid use of persistent variables. And I was already working with my doc writer to clarify this information and add an example of interpolation.

The most up-to-date information on what functions and language features that are supported can be found here. xlsread is not listed, so I didn't think it was supported. Your comments suggest it's supported, but I can't reproduce that on my Mac. To be safe, I suggest using the load command instead. We do document the limitations you're seeing with spline here.

I definitely recommend using a persistent variable rather than a global variable and using load instead of xlsread. I suggesting writing a separate helper function that loads the data using xlsread, creates a piecewise polynomial, and then saves it to the MAT file that you load into the persistent variable.

Finally, if you are able to reproduce the crash, please reach out. That's definitely something we would want to fix ASAP.

[Jim Bosley]

As I mentioned, I was working on using a global variable. global variables seem to be compatible with accelerate/mex. However, after implementing this I found that the ppval function (more specifically the piecewise polynomial, or pp, data structure that ppval uses) aren't compatible with the accelerator function that would speed sbiofit up. So the following pseudo-code works as a function call in a repeat assignment rule in situations not involving accelerate (and mex too, I think), and only calls the xlsread function once:

function glu = glucose_forced(t)
global glupp
if isempty(glupp)
raw = xlsread('my data.xlsx');  % first column is time, second is glucose data
glupp = spline(raw(:,1),raw(:,2));  % or you could use pchip
end
glu = ppval(pp,t);
end

The good news is that the global assignment does seem to be compatible with accelerate (and mex as well, I think). But this highlighted another problem: the pp structure isn't compatible with the accelerate function. The solution to this pp issue that I tried was to use interp1 (1 dimensional interpolation). This looks like:

function glu = forced_glucose(t)
global tgluvec gluvec
if isempty(gluvec)
   raw = xlsread('my data.xlsx');
   tgluvec = raw(:,1);
   gluvec = raw(:,2);
end
glu = interp1(tgluvec,gluvec,t,'spline'): % again, you can specify 'pchip' and probably should.

This does work with a simulate task. That said, I think that pre-processing a piecewise polynomial structure "pp" and using ppval is probably faster, and as I call this function (and several like it) several times each time step, the numerical efficiency may have a slight effect, or more, on run times. I got this second implmentation to work with sbiofit.

The bad news is, after running a script with sbiofit fitting a model using this function to force a parameter variable to change keeps causing MATLAB to crash. As in the dialog box that says wtte "Matlab is stopping now, do you want to sent a message to MATHWORKS about this?" I'll try to get an answer about why this crash is occuring.

To handle the complex QSP models that I'm writing, it is critical that simbiology be able to do things like call simple user-written functions and to have full functionality (that is, the ability to accelerate the model) when those functions are used. Ideally, I'd like accelerate (and mex) to be able to handle 1) persistent variables (is this substantially different than global variables?) and 2) the ability to use functions like spline and pchip to process pp data structures. Perhaps a spline/pchip option to output data in the unmakpp format, and the ability for ppval to use polynomial data in that format? That'd be my suggestion.

If you use a global variable as I do above, make sure that you run

clear global

every time you want to force the function to read the data file.

clear

alone doesn't do it.

Answer 2

[Jim Bosley]

Another update. After I closed Matlab (or rather, it closed itself preemptively!), I opened it up again. I was able to use the model with the interp1 routine (psuedo-code above) with sbiofit, with no error messages regarding accelerate. But when I closed matlab it did crash. Will work with MathWorks folks and report back.

Answer 3

[Jim Bosley]

Arthur, thanks for your note and this information. I'm confused, though, because when I had persistent variables, the error message was thrown in the matlab command line area. I debugged as follows. There was "click here" on the error message, but it didn't work. The error message had a reference to a .mat file in my AppData directory. Loading the .mat file resulted in a matlab variable called report. The report.summary data element referred to an html file. It was actually some sort of matlab format, but when I put it into my browser, Matlab processed it and a very cool error summary popped up, showing me code and line numbers and errors and such. VERY useful. But the error message said that persistent variables weren't supported. And the docs here: https://www.mathworks.com/help/simbio/ug/accelerating-model-simulations-and-analyses.html say "If you have user-defined functions[in Simbiology models], do not use persistent variables in these functions. Persistent variables are not compatible with the functionality used for accelerating simulations." And it goes on to say that if you use persistent variables, you'll get an error message of the form I described above.

After I changed to global variables, I still got an error message, of the same error message format mentioned above. But a different error. It said that the pp struct variable was not consistent with accelerate (and mex? can't recall). So then I changed from using spline/pchip generating pp and using ppeval, to just using interp1 for each call. This (global variables, interp1, and xlsread is the format I use now, and it does not generate errors when I run sbiofit interactively or via script.

Interestingly, xlsread works great in my accelerated functions.

The persistent approach is more satisfying to me, as it's easier, more elegant, and more specific to clear if nothing else (that is, "clear function_name" clears the persistent variables in function_name, but with global variables "clear global" clears everything). Perhaps the error was in the spline/pchip/pp/struct variable area, and persistent variables will work. If I've time I'll try persistent variables with the interp1. That way, if it was spline throwing the error and not persistent variables, we'll know.

I'll check out using load rather than xlsread.

If accelerate really can handle pp struct variables, I'd love to know it: its a more efficient form. But the error messages I got were pretty definitive saying that accelerate wouldn't handle the struct/pp form.

Your help and these exchanges are extremely helpful to me. I'm learning a lot - many thanks. Jim

[Arthur Goldsipe]

Sorry for the confusion... The documentation is incorrect. I can confirm that persistent variables are supported. We only discovered this documentation bug fairly recently, so we have not had time to correct it.

That said, code generation does require you to use all your variables (including persistent variables) more carefully than when writing general MATLAB functions. This was my first time trying to use piecewise polynomials with code generation, so it took me a while to understand the documented constraints. Maybe the easiest way to help you is to give you an example function that is compatible with sbioaccelerate:

    function y = myFun(time)
    persistent pp
    if isempty(pp)
        loaded = load('myData','breaks','coefs');
        pp = mkpp(loaded.breaks, loaded.coefs);
    end
    y = ppval(pp, time);
    end

You just need to use unmkpp to create the variables breaks and coefs from a piecewise polynomial. Then, save those variables to a MAT file. I've attached the one I used for testing this functionality.

Answer 4

[Jim Bosley]

Arthur,

Thanks! This is good news. As I wrote above, the "persistent" approach seems to me to be more elegant than "global".

It is also heartening to know that I can preprocess the forcing data into a pp, and break it out to breaks and coefs, and use those. I'm not sure how much overhead this used in a 35 state model with about 100 reactions, but my coding mentor (John Eaton, a classmate at Texas) taught a pretty intense discipline in coding efficiency.

I use several different data sets with the forcing function, and was just copying the raw (time,datapoint) vectors into the file read by my routines and was letting the simbio routines process the data. This was ok. Now, its clear that it will be more compatible if I preprocess the spline and pchip data to get breaks and coefs. This gives me compatibility AND the speed of preprocessed pp structures. It's a little more manual overhead prior to modeling, but as I understand it, MATLAB has things called scripts I can use to automate the process. ;>)

I'll migrate back to using persistent, and will process my dozen or so datasets into the break/coef mat file format you've suggested. I'll report back on anything I find.

Again, your answers and replies are really helpful.