Make a better API for the BioCro Framework

[leavesandlemmas]

I am making this issue to serve as a reminder of future plans / suggested development. I think we should build a C++ API for biocro.

Here's why:

1. The BioCro developer manual says the BioCro framework should work as a standalone C++ library, so this is an existing goal.
2. Desired changes to the R interface would require making changes to how the framework code works. I think future changes would be facilitated by an API. Desired changes include
   • Making evaluate_module accept a list of modules.
   • We should change how time is specified; allow for saving the state at arbitrary times; remove timestep from the parameters, etc.
3. Building an API might allow us to increase the performance of BioCro as well as increase its functionality. For example:
   • Parameter optimization would be faster if the object BioCro object was instantiated once when defining the objective function, and then persisted for the lifetime of the objective function. Without a way to measure it, I don't know how much this would speed up things (without making the change and testing), but this is the reason using system_derivatives and a R solver is so slow when compared to run_biocro. An API change would also allow interfacing with C++ optimization routines directly.
   • Methods to search module inputs and outputs. I think it would be easier to learn how to use BioCro

Of course, the main reasons against making changes are

1. No one has time. It would take a lot of investment to make happen.
2. Changes would break or change the existing R interfaces which creates additional work.

[justinmcgrath]

It already has an API, which is how R interfaces with C++. This is set up so that the C++ code is independent of R, and it has previously been coupled to OpenSimRoot without using the R interface. That work didn't end up going anywhere though.

I agree that evaluate_module should take a list of modules and we should change how time is specified, although I am not sure the degree to which this requires changing the API. I suspect that the problem here is in the R code, and that the C++ doesn't need to change, but I haven't looked at the details.

Keeping the instance will have minimal effect on performance. Present BioCro is faster than the original version which did everything in C. That indicates that whatever the cost of assembling the object is, it is much, much smaller than running the model. Vaguely considering what gets done during construction, that's probably less work than advancing one time step in the integrator. We solve thousands of time steps, so we'd save a fraction of a percent. There are other reasons to keep the object though, primarily that by separating creating the model and solving it you could use solvers that are not in our C++ library. This would be an API change, but we already interface with the C++ solver directly, which is a constraint, not a goal.

It would be good to have better ways to search search modules and their inputs and outputs. Again, I don't know whether this is an API issue or whether it hasn't been set up in R.

[leavesandlemmas]

agree that evaluate_module should take a list of modules and we should change how time is specified, although I am not sure the degree to which this requires changing the API. I suspect that the problem here is in the R code, and that the C++ doesn't need to change, but I haven't looked at the details.

Here's all of the R /C bindings

extern "C" {
static const R_CallMethodDef callMethods[] = {
    {"R_evaluate_module",                  (DL_FUNC) &R_evaluate_module,                  2},
    {"R_get_all_modules",                  (DL_FUNC) &R_get_all_modules,                  0},
    {"R_get_all_ode_solvers",              (DL_FUNC) &R_get_all_ode_solvers,              0},
    {"R_get_all_quantities",               (DL_FUNC) &R_get_all_quantities,               0},
    {"R_module_creators",                  (DL_FUNC) &R_module_creators,                  1},
    {"R_module_info",                      (DL_FUNC) &R_module_info,                      2},
    {"R_run_biocro",                       (DL_FUNC) &R_run_biocro,                       11},
    {"R_system_derivatives",               (DL_FUNC) &R_system_derivatives,               6},
    {"R_validate_dynamical_system_inputs", (DL_FUNC) &R_validate_dynamical_system_inputs, 6},
    {"R_framework_version",                (DL_FUNC) &R_framework_version,                0},
    {NULL,                                 NULL,                                          0}
};

The only functions which could be used to evaluate modules are R_evaluate_module and R_run_biocro, those both accept external pointers to modules created by R_module_creators. R_evaluate_module accepts a vector of pointers but only evaluates the first one, and it doesn't do any checks. We could have those checks done by validate_dynamical_system_inputs before calling R_evaluate_module (so it would assume that modules form a valid system and are in an order which can be evaluated safely).

Keeping the instance will have minimal effect on performance. Present BioCro is faster than the original version which did everything in C. That indicates that whatever the cost of assembling the object is, it is much, much smaller than running the model.

So why is using system_derivatives with dsolve so much slower then? Even if it's like 100 ms, if run_biocro is called 1000-3000 times during parameter optimization, that works out to be 1-5 min of time spent just constructiong the biocro object even though none of the checks are necessary after the first evaluation. still on the order of minutes of time spent

It would be good to have better ways to search search modules and their inputs and outputs. Again, I don't know whether this is an API issue or whether it hasn't been set up in R.

AFAIK there isn't any way to do this, other than loop through a module library using module_info in R. So it is possible to build something in R that would do that, but if the idea is that BioCro is a standalone C++ library, then I would think it should be possible to ask the biocro object questions like "what module outputs this variable? what module's take it as an input?"

I know that this framework already basically does those calculations. But I don't think it's exposed even in the C++ API.

[justinmcgrath]

The only functions which could be used to evaluate modules are R_evaluate_module and R_run_biocro, those both accept external pointers to modules created by R_module_creators. R_evaluate_module accepts a vector of pointers but only evaluates the first one, and it doesn't do any checks. We could have those checks done by validate_dynamical_system_inputs before calling R_evaluate_module (so it would assume that modules form a valid system and are in an order which can be evaluated safely).

Without looking deeply, this seems like it doesn't need to change the API. The R functions are just not using what is already available. Maybe there are some smallish changes that need to be made though.

So why is using system_derivatives with dsolve so much slower then? Even if it's like 100 ms, if run_biocro is called 1000-3000 times during parameter optimization, that works out to be 1-5 min of time spent just constructiong the biocro object even though none of the checks are necessary after the first evaluation. still on the order of minutes of time spent

I don't know what dsolve is and I've never looked at system_derivatives, so I don't know how it sets up the models. If you're using system_deriviates with something like lsodes, that's not a very fast solver as I recall. It uses a lot of time evaluating things it doesn't need to. Lots of solvers in R are not great, which in part is why we use Boost. In general, I don't see how the fact that system_derivatives takes longer than run_biocro says anything about the cost of constructing the model either. Don't they both construct the model each time they're called? This compares the solvers and whatever cost R has calling the object that system_derivatives returns.

On my computer, running soybean biocro using system.time(with(soybean, run_biocro(initial_values, parameters, soybean_weather[['2005']], direct_modules, differential_modules, ode_solver))) takes 2.85 seconds. If it takes 100 ms to set up the model, that's 3% of the time. That's already slow enough that it's questionable whether that's useful, and 100 ms is a large overestimate. Model construction is assigning pointers and values, operations that take less than a nanosecond, and some simple checks that use cheap operations. Solving the graph is likely the single most time consuming step, and that would take maybe a couple milliseconds. I find it improbable that it's more than 10 milliseconds in total. That's 0.3% of the time. Say you eliminated that 10 milliseconds, running this 1000 times would go from taking 47.5 minutes to 47.33 minutes. Even conceding 100 ms, it would take 45.8 minutes. Even these high side estimates produce differences that are inconsequential, and I am certain that there are better ways to improve performance.

AFAIK there isn't any way to do this, other than loop through a module library using module_info in R. So it is possible to build something in R that would do that, but if the idea is that BioCro is a standalone C++ library, then I would think it should be possible to ask the biocro object questions like "what module outputs this variable? what module's take it as an input?"

This sounds like you could do this in R based on what's available in the API, but it would possibly be useful to instead put the capability directly in the framework and wrap that in R code. Wrapping C++ in R would require passing module libraries as objects though. I'm not entirely sure having this in C++ is useful, since I cannot possibly foresee anyone using this programmatically to build models. It would be restricted to interactive use, so it really only needs to be in R. And putting the logic in R is simpler because it doesn't need to handle passing libraries as external pointers, which adds a lot of complexity. If there's some use case where in C++ on its own someone would want to query what's available in libraries, then I could see spending the work on it. In the absence of that, doing it in R is straightforward and produces something that's useful in most cases for little effort invested.

[leavesandlemmas]

So I have made a version of the framework code and biocro wrapper that allows for a user to call C++ methods of the dynamical system via function calls in R. It's not perfect but it's sufficient to test what the cost of construction of the dynamical_system object is. And @eloch216 ICYMI

The code below is a benchmark, using the simplest ODE we have.

library(BioCro)
library(deSolve)
# creates `dynamical_system` object and returns an external pointer to the object
sys <- 
    create_dynamical_system(
        list(position = 1, velocity = 0),
        list(mass=1, spring_constant=1, timestep=1),
        data.frame(time =seq(0,100,1)),
        list(),
        list("BioCro:harmonic_oscillator")
    )


# this function calls the `calculative_derivative` method on the existing object via the pointer.
f <- function(t, state, param)list(calculative_derivative(sys, 0, state))
run_lsoda_with_handle <- function(){lsoda(c(position=1,velocity=0), seq(0,100,0.5), f)}



# this uses standard biocro; system_derivative returns a function that constructs a `dynamical_system` object at every call of `f2`
f2 <-  system_derivatives(
        list(mass=1, spring_constant=1, timestep=1),
        data.frame(time =seq(0,100,1)),
        list(),
        list("BioCro:harmonic_oscillator"))

run_lsoda_with_system_derivatives <- function(){lsoda(c(position=1,velocity=0), seq(0, 100, 0.5), f2)}

# check output is the same
print(run_lsoda_with_handle())
print(run_lsoda_with_system_derivatives())



# benchmark
microbenchmark::microbenchmark(
    handle = run_lsoda_with_handle(), 
    sys_deriv = run_lsoda_with_system_derivatives(),
    times = 1001)

The benchmark test (on my desktop) produces:

Unit: milliseconds
      expr     min      lq      mean  median      uq     max neval
    handle  2.9222  3.1979  3.476161  3.2826  3.7411 14.1560  1001
 sys_deriv 12.3089 12.7500 13.103679 13.0472 13.3536 25.7183  1001

So the system construction is very slow for this simple example since removing it creates a 3.7x speed up. The soybean model requires more computation per evaluation so the time spent in construction is probably a smaller fraction over total time, but at the same time, the soybean model has many more modules and many direct modules, so the construction process would also require more time.