MCSim under R (Windows)

0.1 Prerequisites

  • R
  • Rtools
  • R Studio (optional)
  • MCSim

This exercise was sourced from the file MCSim under R in GNU MCSim. The testing operation system is Windows 7 with MCSim v5.6.6 and R v3.4.0 installed.

Be sure you’re in the right place. Check .Platform$OS.type == "windows" first!

If you have installed devtools you can use devtools::find_rtools() to check the avaliable of Rtools. Also, we need to chek the Rtools is in the PATH by using

Sys.getenv("PATH")
## [1] "C:\\Program Files\\R\\R-3.5.0\\bin\\x64;C:\\Rtools\\bin;C:\\Windows\\system32;C:\\Windows;C:\\Windows\\System32\\Wbem;C:\\Windows\\System32\\WindowsPowerShell\\v1.0\\;C:\\Program Files (x86)\\Intel\\OpenCL SDK\\2.0\\bin\\x86;C:\\Program Files (x86)\\Intel\\OpenCL SDK\\2.0\\bin\\x64;C:\\Program Files (x86)\\Intel\\Services\\IPT\\;C:\\Program Files\\Intel\\WiFi\\bin\\;C:\\Program Files\\Common Files\\Intel\\WirelessCommon\\;C:\\Program Files (x86)\\Lenovo\\Access Connections\\;C:\\Program Files\\Git\\cmd;C:\\Program Files\\Intel\\WiFi\\bin\\;C:\\Program Files\\Common Files\\Intel\\WirelessCommon\\"

If you can not find the Rtools in the your PATH such as

[1] "c:\\Rtools\\bin;c:\\Rtools\\mingw_64\\bin

You will need to add the Rtools to your path manually, as follow

Sys.setenv(PATH = paste("c:\\Rtools\\bin", Sys.getenv("PATH"), sep=";"))
Sys.setenv(PATH = paste("c:\\Rtools\\mingw_64\\bin", Sys.getenv("PATH"), sep=";"))
Sys.setenv(BINPREF = "c:\\Rtools\\mingw_64\\bin")

Then, use following command to check the function.

Sys.which("gcc")
##                                  gcc 
## "c:\\Rtools\\mingw_64\\bin\\gcc.exe"
system('g++ -v')

If you use git as version control, you can put following files to .gitignore.

*.o
*.exe
*.dll

0.2 Make mod.exe

The following step can compile the simple.model file and create the executable program. The testing model simple.model file looks like this:

#------------------------------------------------------------------------------
# the deSolve example model (simple.model)
#------------------------------------------------------------------------------

States = {y0, y1, y2};

Outputs = {yout};

# Parameters
k1 = 1;
k2 = 1;
k3 = 1;

Dynamics {

  dt(y0) = -k1 * y0 + k2 * y1 * y2; 
  dt(y2) =  k3 * y1 * y1; 
  dt(y1) = -dt(y0) - dt(y2);

  yout = y0 + y1 + y2;

} # End of Dynamics

Events {
  y0 = 1;
}

Roots { # here we need inlining, otherwise 'gout' will not be understood 
  Inline ( gout[0] = y[0] - 0.5; );
}

End.

Make sure that the simple.model and mod folder is put in the working directory!

NOTE: In this case, we don’t need to load libSBML. Therefore, check the config.h file in mod folder. We can see the line 26 in the config.h file like THIS. The #define HAVE_LIBSBML 1 had been disabled. It can prevent the problem when compile to the mod.exe file. (Thanks Yu-Sheng Lin for provding this issue; Mar. 23, 2018).

libSBML = "" # libSBML = "-lsbml" 
# uncomment and execute if you have and want to use libSBML.

# create mod.exe in "mod" folder
system(paste("gcc -o ./mod/mod.exe ./mod/*.c ", libSBML, sep = "")) 

Use model file to check it works.

mName <- "simple.model"

If everything is going well, we can move on to next section.

1 Compile the MCSim model file and solve by deSolve

1.1 Prepare model file

library(deSolve)

# Create .c file and use -R to generate "simple.model_inits.R" initialization file
system(paste("./mod/mod.exe -R ", mName, " ", mName, ".c", sep = "")) 

# compile the MCSim executable, needs to be done each time you modify the model
system (paste0("R CMD SHLIB ", mName, ".c")) # create .o and .dll files
dyn.load(paste(mName, .Platform$dynlib.ext, sep="")) # load simple.model.dll
source ("simple.model_inits.R")

parms_init <- initParms () # define parameter values
Y_init <- initStates (parms) # define initial state values

# set output times, time 0 must be given
times <- c(0, 0.4*10^(0:11))

out.2.1 <- ode(Y_init, times, func = "derivs", parms = parms_init, 
               jacfunc = "jac", dllname = "simple.model", 
               initfunc = "initmod", nout = 1, outnames = "Sum")

plot(out.2.1, log="x",type="b")

We can’t find any change in the dynamic process. It is due to the all initial value equal to 0.

1.2 Customize the parameter and initial values

# Define the parameters and initial values
parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7) 
Y <- c(y1 = 1.0, y2 = 0.0, y3 = 0.0)

# Use the defined parameters and initial values
out.2.2 <- ode(Y, times, func = "derivs", parms = parms,
               dllname = mName, initfunc = "initmod",
               nout = 1, outnames = "Sum")

# Here is the output data frame
head(out.2.2)
##       time        y1           y2         y3 Sum
## [1,] 0e+00 1.0000000 0.000000e+00 0.00000000   1
## [2,] 4e-01 0.9851723 3.386399e-05 0.01479384   1
## [3,] 4e+00 0.9055179 2.240437e-05 0.09445974   1
## [4,] 4e+01 0.7158298 9.185501e-06 0.28416101   1
## [5,] 4e+02 0.4505174 3.222887e-06 0.54947938   1
## [6,] 4e+03 0.1831994 8.942298e-07 0.81679972   1
# Plotting
plot(out.2.2, log="x",type="b")

# Use "which" to select specific variables
plot(out.2.2, log="x",type="b", which = c("y1", "y2", "y3"))

2 Compile the MCSim model file and solve by MCSim

2.1 Prepare in file

In this part, we need to put simple.in file in the working directory. The in file define the initial and parameter values, respectively. Also,we need to set the output time steps and intervals. The source code looks like this:

#------------------------------------------------------------------------------
# the deSolve example model run with GNU MCSim (simple.in)
#------------------------------------------------------------------------------

Integrate (Lsodes, 1E-7, 1E-7, 1);

Simulation {

  y0 = 1;

  # Parameters
  k1 = 0.04;
  k2 = 1e4;
  k3 = 3e7;

  Print (y0, 0, 4e-1, 4e+0, 4e+1, 4e+2, 4e+3, 4e+4, 4e+5, 4e+6, 4e+7, 4e+8, 
         4e+9, 4e+10);
  Print (y1, 0, 4e-1, 4e+0, 4e+1, 4e+2, 4e+3, 4e+4, 4e+5, 4e+6, 4e+7, 4e+8, 
         4e+9, 4e+10);
  Print (y2, 0, 4e-1, 4e+0, 4e+1, 4e+2, 4e+3, 4e+4, 4e+5, 4e+6, 4e+7, 4e+8, 
         4e+9, 4e+10);

  Print (yout, 0, 4e-1, 4e+0, 4e+1, 4e+2, 4e+3, 4e+4, 4e+5, 4e+6, 4e+7, 4e+8, 
         4e+9, 4e+10);

}

End.

The following step compile the simple.model file and create the executable program named mcsim.simple.model.exe.

# set the name of your model
mName = "simple.model"

# Your model executable can be linked with the GNU Scientific Library (GSL, recommended). 
# If you have GSL installed the directives '#define HAVE_LIBGSL 1' and '#define HAVE_LIBGSLCBLAS 1' should appear in the config.h file, in the directory above this one.
# If it does not, put them in, alone at the beginning of a line.  
# If you do NOT have libGSL, those directives should NOT appear in that file. 
# If they do, delete them.

libGSL = "" # libGSL = "-lgsl -lgslcblas" 
# uncomment and execute if you have the GNU Scientific Library


system(paste("./mod/mod.exe ", mName, " ", mName, ".c", sep = "")) 
# Unlike above exercice, you need to remove the -R flag that can only create c file.

Compile the “simple.model.c” to “mcsim.simple.model.exe”. Copy the “sim” file in MCSim folder and put in your working directory. Also copy the config.h file and put in the “sim” folder.

system(paste("gcc -O3 -I.. -I./sim -o mcsim_", mName, ".exe ", mName, ".c ./sim/*.c -lm ", libGSL, sep = ""))

# run the simulation input file
# put its name here
simName = "simple.in"

# run!
system(paste("./mcsim_", mName, ".exe ", simName, sep = ""))

# the output file is called out by default, load it
out.3 <- read.delim("sim.out", skip = 2)

head (out.3)
##    Time       y0          y1        y2 yout
## 1 0e+00 1.000000 0.00000e+00 0.0000000    1
## 2 4e-01 0.985172 3.38640e-05 0.0147937    1
## 3 4e+00 0.905519 2.24048e-05 0.0944587    1
## 4 4e+01 0.715826 9.18549e-06 0.2841650    1
## 5 4e+02 0.450518 3.22290e-06 0.5494780    1
## 6 4e+03 0.183202 8.94234e-07 0.8167970    1
# look
par(mfrow=c(2,2))
plot(out.3[,1], out.3[,2], log="x",type="b", xlab = "time", ylab = "", main = "y1")
plot(out.3[,1], out.3[,3], log="x",type="b", xlab = "time", ylab = "", main = "y2")
plot(out.3[,1], out.3[,4], log="x",type="b", xlab = "time", ylab = "", main = "y3")

Related

comments powered by Disqus