Interact with generated code¶
Interface with the generated code¶
The generated code only represents the mathematical formulation of the model, so in order to find its solution it must be connected to an integration algorithm. This integration is not provided, but some basic functions to interact with the model are. Depending on whether you’re using C++ or Python, the generated code must either be compiled with your solver, or imported as a module.
Build a solver using generated code
This file should be saved as a CMakeLists.txt file for generation with CMake.
cmake_minimum_required(VERSION 3.10.2)
# Delete the previous cache and configured files from this directory before configuring again
# I don't know why it doesn't work ... TODO
file(REMOVE [CMakeCache.txt Makefile cmake_install.cmake modelToSolve.cpp *.h])
file(REMOVE_RECURSE [CMakeFiles])
set(INPUT "" CACHE FILEPATH "Please enter the base name of the generated files to solve (without extension) using the syntax -DINPUT=your_filename")
if("${INPUT}" STREQUAL "")
message(FATAL_ERROR "Please enter the base name of the generated files to solve (without extension) using the syntax -DINPUT=your_filename")
endif()
get_filename_component(INPUTNAME ${INPUT} NAME)
get_filename_component(INPUTDIR ${INPUT} DIRECTORY)
set(PROJECT_NAME solve_${INPUTNAME})
project(${PROJECT_NAME} VERSION 0.1.0)
if(EXISTS "${INPUTDIR}/${INPUTNAME}.c")
configure_file("${INPUTDIR}/${INPUTNAME}.c" "${CMAKE_BINARY_DIR}/modelToSolve.cpp" COPYONLY)
elseif(EXISTS "${INPUTDIR}/${INPUTNAME}.cpp")
configure_file("${INPUTDIR}/${INPUTNAME}.cpp" "${CMAKE_BINARY_DIR}/modelToSolve.cpp" COPYONLY)
endif()
configure_file("${INPUTDIR}/${INPUTNAME}.h" "${CMAKE_BINARY_DIR}/modelToSolve.h" COPYONLY)
configure_file("${INPUTDIR}/${INPUTNAME}.h" "${CMAKE_BINARY_DIR}/${INPUTNAME}.h" COPYONLY)
# TODO This line is a workaround because at present the generator expects the header file to be called "model.h"
# We don't want users to have to modify the generated code, so until there's a better way for users to set
# this in the API, this should stay here. Yes, it's super-clumsy :(
configure_file("${INPUTDIR}/${INPUTNAME}.h" "${CMAKE_BINARY_DIR}/model.h" COPYONLY)
set (PROJECT_SRC
solveGeneratedModel.cpp
modelToSolve.cpp
)
add_executable(${PROJECT_NAME} ${PROJECT_SRC})
message("")
message("1) First use 'make -j' to build the executable")
message("2) Then solve by running: ./${PROJECT_NAME} with the arguments:")
message(" -n step_total")
message(" -dt step_size")
message("")
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 | def module_from_file(input):
# Check the extension is stripped during input.
if input[-3:] != '.py':
module_file = input + ".py"
module_name = input.split("/")[-1]
else:
module_file = input
module_name = ".".join(input.split("/")[-1].split(".")[:-1])
# Import the generated code as a module, and return it.
spec = importlib.util.spec_from_file_location(module_name, module_file)
module = importlib.util.module_from_spec(spec)
sys.modules[module_name] = module
spec.loader.exec_module(module)
return module
|
Full context: solveGeneratedModel.py
Set solution parameters¶
In this very simple example we take the parameters for the integration step size and number of steps from command line arguments. Calling the function below at the start of your program retrieves the command line arguments for use at runtime.
Retrieve command line parameters
std::map<std::string, std::string> processArguments(int argc, char **argv)
{
if (argc == 1) {
std::cout << "Usage:" << std::endl;
std::cout << " -n maxSteps -dt stepSize" << std::endl;
std::cout << " -n the number of steps to take before stopping" << std::endl;
std::cout << " -dt the step size to use" << std::endl;
exit(1);
}
std::map<std::string, std::string> argMap;
std::string value = argv[0];
value.erase(0, 8); // removes the "solve_" from the start of the executable name to get back to the input
argMap["input"] = value;
for (size_t i = 0; i < argc - 1; ++i) {
if (argv[i][0] == '-') {
std::string key = argv[i];
key.erase(0, 1);
value = argv[i + 1];
argMap[key] = value;
i++;
}
}
return argMap;
}
Full context: example_solveGeneratedModel.cpp
def process_arguments(argv):
if (len(argv) == 1):
print("Usage:")
print(" -m generated file to run")
print(" -n the number of steps to take before stopping")
print(" -dt the step size to use")
exit(0)
arg_map = {}
i = 0
while i < len(argv):
if argv[i][0] == '-':
key = argv[i][1:]
value = argv[i + 1]
arg_map[key] = value
i += 1
else:
i += 1
# Cleaning up the inputs to save in the right form
error_string = ''
try:
arg_map['m'][-3:] == ".py"
arg_map['m'] = arg_map['m'][:-3]
except:
error_string += "/n - missing argument: -m file to run"
try:
arg_map['n'] = int(arg_map['n'])
except:
error_string += "/n - missing argument: -n number of steps to take"
try:
arg_map['dt'] = float(arg_map['dt'])
except:
error_string += "/n - missing argument: -dt step size"
if error_string != "":
print(error_string)
exit(1)
return arg_map
Full context: example_solveGeneratedModel.py
The specific parameters related to the model are available within the generated code, and can be retrieved from their structures as shown below.
Retrieve model parameters
// Inside the generated code are structures with information about the
// model and its dimensions. These are:
// - VOI_INFO: a dict with the name, units, and component of the variable of integration,
// - STATE_INFO: a list of dicts for the state variables,
// - VARIABLE_INFO: a list of dicts for the non-state variables.
std::cout << " VARIABLE OF INTEGRATION (units) " << std::endl;
std::cout << " " << VOI_INFO.name << " (" << VOI_INFO.units << ")" <<std::endl
<< std::endl;
std::cout << " STATE VARIABLES (units) " << std::endl;
std::cout << "-------------------------------------------------------------" << std::endl;
for (size_t i = 0; i < STATE_COUNT; ++i) {
std::cout << " " << STATE_INFO[i].name << " (" << STATE_INFO[i].units << ")" << std::endl;
}
std::cout << std::endl;
Full context: example_solveGeneratedModel.cpp
# Retrieve model module from the generated code file.
model = module_from_file(args['m'])
# Inside the 'model' module are structures with information about the
# model and its dimensions. These are:
# - VOI_INFO: a dict with the name, units, and component of the variable of integration,
# - STATE_INFO: a list of similar dicts for the state variables,
# - VARIABLE_INFO: a list of similar dicts for the non-state variables.
print(' VARIABLE OF INTEGRATION (units)')
print('--------------------------------------------------------------------')
print(' {} ({}, {})'.format(model.VOI_INFO['name'],
model.VOI_INFO['units'],
dt))
print(' {} steps'.format(n))
print()
Full context: example_solveGeneratedModel.py
Initialise solution variables¶
As well as setting the number and size of the steps, the basic integration algorithm needs to store the state variables and their rates, as well as any other variables that are independent of the integration.
Allocate solution arrays
// Call module functions to construct and initialise the variable arrays.
// Note that both the rates and the states arrays have the same dimensions,
// so it's possible to call the createStatesArray() function for both.
double time = 0.0;
auto myStateVariables = createStatesArray();
auto myRates = createStatesArray();
auto myVariables = createVariablesArray();
Full context: example_solveGeneratedModel.cpp
# Call module functions to construct the variable arrays.
# Note that both the rates and the states arrays have the same dimensions,
# so it's possible to call the create_states_array() function for both.
time = 0.0
my_variables = model.create_variables_array()
my_state_variables = model.create_states_array()
my_rates = model.create_states_array()
Full context: example_solveGeneratedModel.py
Note that since the number of state variables must be identical to the number of rates, the same create states array function can be used to allocate storage for both of them.
Following allocation, the states and constants can be initialised using the function provided.
Compute the starting values of variables¶
There are two kinds of variables stored in the generated code: those that require updating each step (the rates, but also any other variables which these depend upon); and those which do not affect the progress of the solution (everything else). The former kind is updated when the compute rates function is called, and the latter when the compute variables function is called. Constants clearly need be calculated only once, using the compute computed constants function.
In this example we compute everything at the initial point so as to be able to output the starting values to the terminal.
Initialise variables
// Make use of the access functions provided to initialise the variable arrays.
initialiseStatesAndConstants(myStateVariables, myVariables);
computeComputedConstants(myVariables);
computeRates(time, myStateVariables, myRates, myVariables);
computeVariables(time, myStateVariables, myRates, myVariables);
Full context: example_solveGeneratedModel.cpp
# Compute the parameters which require it, including the rates and variable values.
model.initialise_states_and_constants(my_state_variables, my_variables)
model.compute_computed_constants(my_variables)
model.compute_rates(0, my_state_variables, my_rates, my_variables)
model.compute_variables(0, my_state_variables, my_rates, my_variables)
print(' STATE VARIABLES (units, initial value)')
print('--------------------------------------------------------------------')
for i in range(0, model.STATE_COUNT):
print(' {} ({}, {})'.format(model.STATE_INFO[i]['name'],
model.STATE_INFO[i]['units'],
my_state_variables[i]))
print()
print(' VARIABLES (units, initial value)')
print('--------------------------------------------------------------------')
for v in range(0, model.VARIABLE_COUNT):
print(' {} ({}, {})'.format(model.VARIABLE_INFO[v]['name'],
model.VARIABLE_INFO[v]['units'],
my_variables[v]))
Full context: example_solveGeneratedModel.py
Prepare files for output¶
The solution is to be written directly to an output file during the iteration process, and this step is simply the preparation for that. The file is created and opened, and the columns labelled with information from the model.
Prepare for output
// Prepare a file for writing during the solution process.
std::cout << " INITIAL CONDITIONS" << std::endl;
std::cout << "-------------------------------------------------------------" << std::endl;
for (size_t i = 0; i < STATE_COUNT; ++i) {
std::cout << " " << STATE_INFO[i].name << "(" << VOI_INFO.name << " = 0) = " << myStateVariables[i] << std::endl;
}
for (size_t i = 0; i < VARIABLE_COUNT; ++i) {
std::cout << " " << VARIABLE_INFO[i].name << "(" << VOI_INFO.name << " = 0) = " << myVariables[i] << std::endl;
}
std::cout << std::endl;
std::string outFileName = args["input"] + "_solution.txt";
std::ofstream outFile(outFileName);
// Header line for output file
outFile << "iteration";
outFile << "\t" << VOI_INFO.name << " (" << VOI_INFO.units << ")";
for (size_t s = 0; s < STATE_COUNT; ++s) {
outFile << "\t" << STATE_INFO[s].name;
}
for (size_t s = 0; s < VARIABLE_COUNT; ++s) {
outFile << "\t" << VARIABLE_INFO[s].name;
}
outFile << std::endl;
// Initial conditions in output file
outFile << 0;
outFile << "\t" << 0;
for (size_t s = 0; s < STATE_COUNT; ++s) {
outFile << "\t" << myStateVariables[s];
}
for (size_t s = 0; s < VARIABLE_COUNT; ++s) {
outFile << "\t" << myVariables[s];
}
outFile << std::endl;
Full context: example_solveGeneratedModel.cpp
# Prepare to write output to a file during the solution process.
row = 'iteration\t{}({})'.format(
model.VOI_INFO['name'], model.VOI_INFO['units'])
for s in range(0, model.STATE_COUNT):
row += '\t{}({})'.format(model.STATE_INFO[s]
['name'], model.STATE_INFO[s]['units'])
for s in range(0, model.VARIABLE_COUNT):
row += '\t{}({})'.format(model.VARIABLE_INFO[s]
['name'], model.VARIABLE_INFO[s]['units'])
row += '\n'
write_file_name = '{}_solution.txt'.format(args['m'])
write_file = open(write_file_name, 'w')
write_file.write(row)
Full context: example_solveGeneratedModel.py
Perform the integration¶
There are myriad stepping schemes for numerical integration, but here we use the very basic Euler method. At each step, new rates are calculated by calling the compute rates function. The states are then extrapolated from the rates using the Euler approximation step. In this example we are writing the output at every step, so we also need to update the variables too using the compute variables function, and everything is written to the output file.
Perform the integration and output
// Numerically integrate the state variables using the Euler method to step through the solution.
// Solution columns in output file
for (size_t step = 1; step < stepCount; ++step) {
time = step * stepSize;
computeRates(time, myStateVariables, myRates, myVariables);
outFile << step << "\t " << time;
for (size_t s = 0; s < STATE_COUNT; ++s) {
myStateVariables[s] = myStateVariables[s] + myRates[s] * stepSize;
outFile << "\t" << myStateVariables[s];
}
// The variables in the "myVariables" array are those which do not affect the calculation
// of rates or state variables. They only need to be computed when outputting the
// results of a timestep: if you're not saving every timestep, then you can skip this
// until you are.
computeVariables(time, myStateVariables, myRates, myVariables);
for (size_t s = 0; s < VARIABLE_COUNT; ++s) {
outFile << "\t" << myVariables[s];
}
outFile << "\n";
}
outFile.close();
// Housekeeping.
deleteArray(myStateVariables);
deleteArray(myVariables);
deleteArray(myRates);
Full context: example_solveGeneratedModel.cpp
# Numerically integrate using Euler steps to solve the model.
for step in range(0, n):
time = step * dt
model.compute_rates(time, my_state_variables, my_rates, my_variables)
# Formatting for output.
row = '{}\t{}'.format(step, time)
for s in range(0, model.STATE_COUNT):
my_state_variables[s] = my_state_variables[s] + \
my_rates[s] * dt
row += '\t{}'.format(my_state_variables[s])
# Note that the variables in the my_variables array are those which
# are independent of the integration: thus, they only need to be
# computed at timesteps where the solution is to be written to the
# output. For large simulations, this may not be every integration
# timestep.
model.compute_variables(
time, my_state_variables, my_rates, my_variables)
# Output the solution.
for s in range(0, model.VARIABLE_COUNT):
row += '\t{}'.format(my_variables[s])
row += '\n'
write_file.write(row)
write_file.close()
Full context: example_solveGeneratedModel.py