Last week I posted a C++ implementation of basic cancer stem cell driven tumor growth model. About 100 lines of a code allowed to perform simulation, but there was nothing about exporting the results, doing more complicated analysis, visualization, or performing data fitting. Of course each of those tasks can be performed by writing another large piece of the code.

Two weeks ago I posted a MATLAB code for the very similar model, which was quite fast but a lot slower than C++ implementation. On the other hand it was straightforward to visualize the tumor and it’s not much of a hassle to perform more complicated analysis or data fitting using variety of MATLABs built-in functions.

Each of those two approaches has their own advantages: C++ is fast, MATLAB has a lot of built-in in tools.

Today I will show how to take the best out of both approaches, i.e. I’ll show how to easily wrap our C++ code into MATLAB. The whole C++ based simulation will be invoked from MATLAB as a function and results will be visible straight from it.

First we need to add necessary libraries.

#include "mex.h" #include <matrix.h>

Then we need to modify our original definition of global variables (see previous post) so they can be adjusted based on the input parameters (for example lattice has to be dynamically allocated).

int N; bool* lattice; vector<cell> cells; char pmax; double pDiv, alpha, ps, pmig; int* indcNeigh;

Then the only thing that we need to do is to delete old main() function and use mexFunction() instead. In the mexFunction() we need to read first the parameters that will be passed from MATLAB (input parameters to the function). Then we need to assign the memory for the lattice and set other auxiliary variables. Finally, after performing the simulation we need to write variables to the output, so the MATLAB can see the results.

void mexFunction( int nlhs, //number of output arguments mxArray *plhs[], //output arrays int nrhs, //number of input arguments const mxArray *prhs[] //input arrays ) { /*READING INPUT PARAMETERS*/ double *par=mxGetPr(prhs[0]); //model parameters, first input variable double *settings=mxGetPr(prhs[1]); //simulation settings, second input variable int Nsteps = settings[0]; //number of simulation steps to perform N=(int)settings[1]; //size of the lattice pmax = par[0]; //divisions before proliferative capacity exhaustion pDiv = par[1]; //division probability alpha = par[2]; //spontaneous death probability ps = par[3]; //probability of symmetric division pmig = par[4]; //probability of migration /*CREATING LATTICE AND AUXILARY VARIABLES*/ //creating lattice with a given size lattice = new bool [N*N]; fill_n(lattice, N*N, false); indcNeigh = new int [8];//cell neighborhood //assigning values indcNeigh[0] = -N-1; indcNeigh[1] = -N; indcNeigh[2] = -N+1; indcNeigh[3] = -1; indcNeigh[4] = 1; indcNeigh[5] = N-1; indcNeigh[6] = N; indcNeigh[7] = N+1; /*OLD PIECE OF MAIN() FUNCTION*/ srand(time(NULL)); //initialize random number generator initialize(); //initialize CA simulate(Nsteps); //simulate CA /*WRITING TO OUTPUT VARIABLES*/ //1) writing the lattice as a logical matrix plhs[0] = mxCreateLogicalMatrix(N,N); //lattice mxLogical *ptr = mxGetLogicals(plhs[0]); for(int i=0; i<N*N; i++) ptr[i] = lattice[i]; //2) writing the cells vector as a vector plhs[1] = mxCreateDoubleMatrix(1,cells.size()*3,mxREAL); double* ptr2 = mxGetPr(plhs[1]); for(int i=0; i<cells.size(); i++) { ptr2[3*i] = cells.at(i).place; ptr2[3*i+1] = cells.at(i).p; ptr2[3*i+2] = cells.at(i).is_stem; } }

After making all of the changes we can compile our .cpp file using the MATLAB command “mex CA.cpp”. Remember all other pieces of C++ associated with model itself remain the same – there is no need to change anything.

From now on we can perform our simulation in MATLAB environment by invoking the following command

[lattice, cells] = CA([pmax, pdiv, alpha, ps,pmig],[nSteps, N]);

and further use of all MATLABs benefits, especially data fitting features. At the same time all of the simulations will be performed at C++ speed! Imagine using Parallel Toolbox in MATLAB to harness multiple CPUs with only few lines of the code.

[…] weeks ago I’ve shown here how to wrap your C++ in order to use it from within MATLAB as a function without loosing the C++ […]

LikeLike

[…] Of course, the trade-off is that it won’t be that easy to read our simulation files to other programs. However, we can always write wrapper functions that will allow other programs to interpret our files (see e.g. this post). […]

LikeLike