specification.qmd

```{r,echo=FALSE,message=FALSE}
source("setup.R")
```

# Model specification {#sec-model-specification}

This chapter details the mrgsolve model specification format.  

## How / where to write a model

There are two ways to write your model:

1. Code in a __separate file__ and source into your R script 
2. Code __inline__ as a character string already in your R script

We recommend method `1` (separate file) for any non-trivial modeling work. 
Method `2` is handy for quickly coding a model and you'll also see us using 
that approach frequently when demonstrating how to use mrgsolve.

### Separate file {#sec-spec-separate-file}

For most applications, you will want to put your model code in a standalone
file. This file can have any extension, but there is some special behavior 
when you use the `.cpp` extension. We'll show you both here. 

#### Using `.cpp` file extension

Open a text editor and type the model code into a file with name that has the
format `<model-name>.cpp`.  This filename format identifies a "name" for your
model (`<model-name>`, the "stem" of the file name).  Note: this whole file will
be read and parsed, so everything in it must be valid mrgsolve model
specification elements.  

Use the `mread()` function to read and parse this file.  For the model 
called `mymodel` saved in `mymodel.cpp` (in the current working directory), 
issue the command:

```{r,eval=FALSE}
mod <- mread("mymodel")
```

`mread()` returns a model object from which you can simulate.

#### Using any file extension

You can use any file extension to save your mrgsolve model. Rather 
than using `.cpp`, you might use `.mod`. When you are using an extension
other than `.cpp`, pass in the entire file name, both stem and 
extension. If you have the model code saved in `mymodel.mod` then 
you'd call

```{r, eval = FALSE}
mod <- mread("mymodel.mod")
```

to load that model. 

#### The model project directory

The second argument to `mread()` is `project`. This is the directory where
mrgsolve will look for your model file (passed as the first argument to
`mread()`). 

If your `myproject.mod` file is located in the `models` directory, then 
you could load that file with 

```{r, eval = FALSE}
mod <- mread("mymodel.mod", project = "models")
```

Alternatively, you can just past the entire path to the model file

```{r, eval = FALSE}
mod <- mread("models/mymodel.mod")
```


### Inline / `code` \ {#sec-spec-inline}
Often it is more convenient to write a model right in your `R` script.  The 
model might look something like this:

```{r}
code <- '
$PARAM CL = 1, VC = 20
$PKMODEL ncmt=1
'
```

Here, we created a character vector of length 1 and saved it to the `R` object
called `code`.  The name of this object is irrelevant.  But `code` will be 
passed into mrgsolve as the model definition.  When mrgsolve gets a model 
like this along with a "name" for the model, mrgsolve will write the code to a 
file called `<model-name>.cpp` and read it right back in as if you had typed 
the code into this file (@sec-spec-separate-file). 

To parse and load this model, use the `mcode()` command:

```{r,eval=FALSE}
mod <- mcode("mymodel", code)
```

`mcode()` is a convenience wrapper for `mread()`.  `mcode` writes the code to
`mymodel.cpp` in `tempdir()`, reads it back in, compiles and loads.

The `mcode` call is equivalent to:
```{r, eval=FALSE}
mod <- mread("mymodel", tempdir(), code)
```

For help, see `?mread` , `?mcode` in the `R` help system after loading 
mrgsolve.

### Comments

You can comment out code in the model specification file with two 
front slashes `//`

```{c, eval = FALSE}
$MAIN

double a = 1.234;
// double b = 5.678; 
```

This comment sequence is what is used in `C++` code, but you can use it 
_anywhere_ in the mrgsolve model specification file.

## Code blocks {#sec-code-blocks}

### About code blocks


__Block identifier__

Different types of code are organized in the model specification file and 
separated by block identifiers. There are two ways to formulate block 
identifiers that can be used in mrgsolve.  In the first type, a dollar-sign is 
placed at the start of the block name

```{c,eval=FALSE}
$BLOCKNAME
<block-code>
```

For example, a block of parameters would be

```{c, eval=FALSE}
$PARAM
CL = 1
```

The second way to write this is with brackets

```{c, eval=FALSE}
[ BLOCKNAME ] 
<block-code>
```

There is no functional difference between the dollar-sign notation and the 
brackets.  When model specification code is saved into a file with a `.cpp` 
extension, the code editor may make certain assumptions about formatting or
styling the code.  Using brackets will most-likely work better with the 
editor in that case.  

Block identifiers are case-insensitive so all of these also work

```{c, eval = FALSE} 
$param 
CL = 1
```

```{c, eval = FALSE}
[ param ] 
CL = 1
```

Users are free to include block code on the same line as the block identifier,
but must include a space after the identifier. For example, the parser will
recognize `$PARAM CL = 1` but not `$PARAMCL=1` as parameters.


__Block syntax__
Different blocks may require different syntax.  For example, code written in
`$PARAM` will be parsed by the `R` parser and will generally need to adhere to
`R` syntax requirements.  On the other hand, code in `$MAIN`, `$ODE`, and
`$TABLE` will be used to make functions in `C++` and therefore will need to be
valid `C++` code, including terminal `;` on each line.


__Block options__
Options may be specified on some code blocks that signal how the code is 
to be parsed or used in the simulation.  

Block options can be boolean in which case you indicate the option name 
with nothing following.  For example, to indicate `block = TRUE` for an 
OMEGA matrix, write

```{c, eval = FALSE}
$OMEGA @block 
1 2 3
```

Here `@block` indicates `block = TRUE`. Starting with mrgsolve 1.0.0, block 
options can be negated by writing `!` between `@` and the option name. To 
request a non-block (diagonal) OMEGA matrix

```{c, eval = FALSE}
$OMEGA @!block 
1 2 3
```

But note well that diagonal is the default configuration for OMEGA, so it is 
never necessary to write this. 

Other non-boolean block options state the name and then a value to be assigned
to that name. For example, if we wanted to assign labels to the ETAs from 
an OMEGA block, we'd write

```{c, eval = FALSE}
$OMEGA @labels ECL EVC
1 2
```

This essentially sets labels equal to the vector `c("ECL", "EVC")`.

We note two specific boolean options (`@object` and `@as_object`) that have 
similar function across multiple blocks.  

### Programmatic or bulk initialization 

The following describes syntax for initializing blocks as R objects using
R code. This can be helpful, say, when you need to initialize a 50x50
OMEGA matrix or a series of systematically named parameters. We will describe
the `@object` and `@as_object` options that are available on select blocks.

The `@object` option lets you name an object defined in `$ENV` to use to 
instantiate the block data.  For example, to specify a series of parameters
using the `@object` option, you'd write

```{c, eval = FALSE}
$ENV
params <- list(CL = 1, V = 20, KA = 1.2)

$PARAM @object params
```

This tells mrgsolve that there is an object called `params` that it is in 
the `$ENV` environment and mrgsolve will use that to define the names and 
values. This is a trivial example to show how a simple series of parameters
could be defined. However, the intended use for this functionality is to allow
efficient creation of a large, systematically-named series of parameters in 
a large model. 

The `@as_object` option is a boolean option that tells mrgsolve that the 
block code will actually return the object (rather than asking mrgsolve to 
look in `$ENV` for the object). The equivalent specification for the block 
above would be:

```{c, eval = FALSE} 
$PARAM @as_object
list(CL = 1, V = 20, KA = 1.2)
```

The following blocks contain both the `@object` and `@as_object` options:

- `$PARAM`
- `$INPUT`
- `$THETA`
- `$CMT`
- `$INT`
- `$OMEGA`
- `$SIGMA`

Please see the specific block documentation for more details on the specific
type of object that should be returned when this syntax is invoked. 


### `$PROB` \ {#sec-block-prob}

__Syntax__: text   
__Multiple allowed__: yes   
__Options__: `@annotated`, `@covariates`

Use this block to make notes about the model.  There are no restrictions on the 
text that gets entered here.  mrgsolve does not routinely process the text 
in any way, except when rendering the model as a document.  Frequently, we 
write the text here in markdown format so that it will render nicely in the
model document.  But this is completely optional. 

See the annotated model in @sec-topic-annotated for an example. 

### `$PARAM` \ {#sec-block-param}

__Syntax__: R   
__Multiple allowed__: yes   
__Options__: `@annotated`, `@covariates`, `@tag`, `@input`, `@object`, `@as_object`

Define the parameter list in the current model. Parameters are names associated
with values that can be used throughout the model.  A value must be given 
for every parameter name.  Names (and numbers) of parameters must be set at the
time the model is compiled, but parameter values may be updated without 
re-compiling the model. 

The `@covariates` option allows you to tag a collection of parameters as 
"covariates".  It does not change the functionality of the model or simulation
workflow in any way, but allows you to get that list of covariate names out
of the model object.

The `@input` option adds a tag to the parameters in the block so that input
data sets will be checked for these parameters when `check_data_names()` is
called. Use the `@tag` option and supply a user-defined tag that can also be 
used for checking data sets and parameters. 

Example:
```{c, eval=FALSE}
[ PARAM ] CL = 1, VC = 20, KA = 1.2
KM = 25, VMAX = 400, FLAG = 1, WT = 80
SEX = 0, N = sqrt(25)
```

Note that the "values" in the parameter list will get evaluated by the R 
interpreter and the evaluation will happen inside the environment which is 
defined by the `$ENV` block. For example

```{c, eval = FALSE}
$ENV MWT = 1.2

$PARAM 
WT = 80 * 2.2
MASS = 600 * MWT
```

Annotated example:
```{c, eval=FALSE}
[ PARAM ] @annotated
CL :   1 : Clearance (L/hr)
VC :  20 : Volume of distribution (L)
KA:  1.2 : Absorption rate constant (1/hr)
```

See the `popex` model for an example of using `@covariates`

```{r, include = FALSE, message = FALSE}
obj <- as.list(mod <- modlib("popex", compile = FALSE))
```

```{r, echo = FALSE}
blocks(mod, PARAM)
```

then

```{r, eval  = FALSE}
mod <- modlib("popex")
```

```{r}
as.list(mod)$covariates
```

Notes:

* Multiple blocks are allowed
* Values are evaluated by the `R` interpreter

The `@object` and `@as_object` options should name or return a named list of 
parameters.

See also:  @sec-block-input, @sec-block-theta and @sec-block-fixed.

See `?param` in the `R` help system after loading mrgsolve.


### `$INPUT` \ {#sec-block-input}

__Syntax__: R   
__Multiple allowed__: yes   
__Options__: `@annotated`, `@covariates`, `@tag`, `@object`, `@as_object`

This block operates just like $PARAM (@sec-block-param), in that it is a way
to introduce parameters in to the model. Additionally, the `@input` tag is 
added to the model parameters and input data sets will be checked for these 
parameters when `check_data_names()` is called. 

See also:  @sec-block-param, @sec-block-theta,  @sec-block-fixed.


### `$FIXED`  \ {#sec-block-fixed}

__Syntax__: R   
__Multiple allowed__: yes   
__Options__: `@annotated`, `@object`, `@as_object`


Like `$PARAM`, `$FIXED` is used to specify `name=value` pairs.  Unlike 
`$PARAM`, however, the values associated with names in `$FIXED` are not able 
to be updated.  

By default, names in `$FIXED` are associated with their value through a 
`C++` preprocessor `#define` statement. 

Usually, `$FIXED` is only used when there are a very large number of 
parameters ($>$ 100 or 200).  When some of these parameters never need 
to be updated, you can move them to a `$FIXED` block to get a modest gain 
in efficiency of the simulation.

Items in `$FIXED` will not be shown when parameters are queried. 

Example:
```{c,eval=FALSE}
[ PARAM ] CL = 2, VC = 20

[ FIXED ]
g = 9.8
```


Annotated example:
```{c, eval=FALSE}
$FIXED @annotated
g : 9.8 : Acceleration due to gravity (m/s^2)
```

See also: @sec-block-param and @sec-block-theta.

Notes:

* Multiple blocks are allowed
* Values are evaluated by the `R` interpreter


### `$CMT`  and  `$INIT`

__Syntax__: text   
__Multiple allowed__: yes   
__Options__: `@annotated`, `@object`, `@as_object`

Declare the names of all compartments in the model. 

* For `$CMT` give the names of compartments; initial values are assumed to be 0
* For `$INIT` give the name and initial value for all compartments

Note that both `$CMT` and `$INIT` declare compartments, so any compartment name
should get declared in either `$CMT` or `$INIT`, but never both.

Examples:
```{c,eval=FALSE}
[ CMT ] GUT CENT RESPONSE
```

```{c,eval=FALSE}
[ INIT ] GUT  = 0, CENT = 0, RESPONSE = 25
```

Annotated examples:
```{c, eval=FALSE}
[ CMT ] @annotated
GUT      : Dosing compartment (mg)
CENT     : Central PK compartment (mg)
RESPONSE : Response
```

```{c, eval=FALSE}
$INIT @annotated
GUT      :   0 : Dosing compartment (mg)
CENT     :   0 : Central PK compartment (mg)
RESPONSE :  25 : Response
```

The `@object` and `@as_object` options should name or return a named list of 
compartments and initial values when used in `$INIT` and a character vector
of compartment names when used in `$CMT`.

See `?init` in the `R` help system after loading mrgsolve.

### `$MAIN`  \ {#sec-block-main}

__Syntax__: C++   
__Multiple allowed__: no   
__Options__: `@check_modeled_infusions`

This code block has two main purposes:

* Derive new algebraic relationships between parameters, random, effects and 
  other derived variables
* Set the initial conditions for model compartments

For users who are familiar with NONMEM, `$MAIN` is similar to `$PK`.

`$MAIN` is wrapped into a `C++` function and compiled / loaded by mrgsolve.  
The `MAIN` function gets called just prior to advancing the system from the 
current time to the next time for each record in the data set.  `$MAIN` also 
gets called several times before starting the problem (`NEWIND == 0`) and 
just prior to simulating each individual (`NEWIND == 1`).  Finally, `$MAIN` 
gets called every time the model initial conditions are queried with `init()`.  

New variables may be declared in `$MAIN`.  See @sec-new-variables for details.

Examples:

```{c,eval=FALSE}
[ CMT ] CENT RESP

[ PARAM ] KIN = 100, KOUT = 2, CL = 1, VC = 20

[ MAIN ]

RESP_0 = KIN/KOUT;

double ke = CL/VC;

```


Note that `@check_modeled_infusions` block option defaults to true. To turn off
warnings from mrgsolve when checking modeled infusions against the data set 
`RATE`, use the negation: `@!check_modeled_infusions`. mrgsolve will always 
error when an infusion resolves to a negative rate or duration.


### `$PK` \ {#sec-block-pk}

This is an alias for `$MAIN`.


### `$ODE` \ {#sec-block-ode}

__Syntax__: C++   
__Multiple allowed__: yes
__Options__: `@param`

Use `$ODE` to define model differential equations.  For all compartments assign 
the value of the differential equation to `dxdt_CMT` where `CMT` is the name 
of the compartment.  The `dxdt_` equation may be a function of model 
parameters (via `$PARAM`), the current value of any compartment (`CMT`) 
or any user-derived variable. 

For example:
```{c,eval=FALSE}
[ CMT ] GUT CENT

[ ODE ]
dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - KE*CENT;

```

It is important to make sure that there is a `dxdt_` expression defined for 
every compartment listed in `$CMT` or `$INIT`, even if it is `dxdt_CMT = 0;`

The `$ODE` function is called repeatedly during a simulation run.  So it is 
wise to do as many calculations as possible outside of `$ODE`, usually in 
`$MAIN`.  But remember that any calculation that depends on an amount in a
compartment and helps determine the `dxdt_` expression in a model must be 
written in `$ODE`.


New variables may be declared in `$ODE`.  See @sec-new-variables
for details.

For example:
```{c,eval=FALSE}
$CMT CENT RESP
$PARAM VC = 100, KE = 0.2, KOUT = 2, KIN = 100
$ODE
double CP = CENT/VC;
double INH = CP/(IMAX+CP)

dxdt_CENT = -KE*CENT;
dxdt_RESP =  KIN*(1 - INH) - RESP*KOUT;

```

If the model needs to refer to the current time, use the `SOLVERTIME` variable.

Notes:

* `$ODE` is written in `C++` syntax; every line must end in `;`
* There may be only one `$ODE` block in a model

### `$DES` \ {#sec-block-des}

This is an alias for `$ODE`.


### `$TABLE` \ {#sec-block-table}

__Syntax__: C++   
__Multiple allowed__: yes 

Use `$TABLE` to interact with parameters, compartment values, and other
user-defined variables __after__ the system advances to the next time.

For example:
```{c,eval=FALSE}
[ TABLE ]
double CP = CENT/VC;
```

__NOTE__
mrgsolve formerly had a `table()` macro for inserting derived values into
simulated output. This macro has been deprecated.  The only way to insert
derived values into the simulated output is via `$CAPTURE`.

__NOTE__
When variables are marked for capture (see @sec-block-capture), the values of
those variables are saved at the __end__ of the `$TABLE` function.  This process
is carried out automatically by mrgsolve and therefore requires no user
intervention.

### `$ERROR` \ {#sec-block-error}

This is an alias for `$TABLE`.

### `$EVENT` \ {#sec-block-event}

__Syntax__: C++     
__Multiple allowed__: no      
__mrgsolve version__: >= 1.5.2   

Use the `$EVENT` block to write code to implement doses or other modeled events
inside your model file. `$EVENT` functions just like `$ERROR` (or `$TABLE`),
but it gets called just prior to `$ERROR`. This lets you set up modeled doses
or events for execution "now" and also seem them in the simulated output on 
that same record. For example


```{r}
#| eval: false
$EVENT 
if(NEWIND > 1) return;

evt::ev dose = evt::bolus(250, 2);
evt::ii(dose, 12);
evt::ss(dose, 1); 
evt::addl(dose, 9)

self.push(dose);
```

**A subtle difference between `$ERROR` and `$EVENT`**

Coding modeled events in `$EVENT` rather than `$TABLE` will give a very limited
and subtle difference in output that you will most likely see when you
administer a bolus dose to a compartment you are watching (e.g. a bolus to the
central compartment). Using `$EVENT` will let you see the result of the dose
immediately - on the record where the dose code is executed; when `$TABLE` is
used, you will only see the result of that dose on subsequent output records.
You'd see a difference between `$TABLE` and `$EVENT` behavior when you bolus
into a depot compartment too, but you probably aren't noticing because you are
likely monitoring the central compartment, not the depot. Regardless, the
`$TABLE` / `$EVENT` difference is just bookkeeping for the current record; you
might see that bookkeeping difference on the current record, but the difference
will not factor into how the system advances beyond that record.

**The key point for using `$EVENT`**

For the time being, the important thing to remember when you want to dose
"now" (not some time in the future) is those interventions should be coded into 
`$EVENT` (or `$ERROR`) and they will happen right after the system advances 
in time on the current record, not prior to advancing in time for the next 
record. I recommend you always use `$EVENT` or `$TABLE` to execute doses or 
other events from inside your model. This distinction only applies to doses
happening "now"; there is no difference in behavior when the dose is scheduled
to happen at a future time.


### `$PREAMBLE` \ {#sec-block-preamble}

__Syntax__: C++   
__Multiple allowed__: no   

This is the fourth C++ code block.  It is called once in two different settings
(@sec-section-sequence):

1. Immediately prior to starting the simulation run
1. Immediately prior to calling `$MAIN` when calculating initial conditions


`$PREAMBLE` is a function that allows you to set up your C++ environment.  It is
only called one time during the simulation run (right at the start).  The code
in this block is typically used to configure or initialize C++ variables or data
structures that were declared in `$GLOBAL`.

For example:

```{c,eval=FALSE}
[ PLUGIN ] Rcpp

[ GLOBAL ]
namespace{
  Rcpp::NumericVector x;
}

[ PREAMBLE ]
x.push_back(1);
x.push_back(2);
x.push_back(3);

[ MAIN ]
<some code that uses x vector>
```

In this example, we want to use a numeric vector `x` and declare it in `$GLOBAL`
so that we can use it anywhere else in the code (the declaration is also made in
an unnamed namespace to ensure that the variable is local to the model file).
Then, in `$PREAMBLE`, we put 3 numbers into the vector and we use `x` in
`$MAIN`.  Since `$MAIN`, `$TABLE` and (especially) `$ODE` are called repeatedly
as the simulation run progresses, we put the initialization of `x` in `$PREAMBLE`
to make sure the initialization of `x` only happens once.

Notes:

* `$PREAMBLE` is written in `C++` syntax; every line must end in `;`
* There may be only one `$PREAMBLE` block in a model
* Like `$MAIN`, `$ODE` and `$TABLE`, `double`, `int` and `bool` variables 
  initialized in `$PREAMBLE` are actually initialized for global (within the 
  model file)

See also: @sec-section-sequence and @sec-block-plugin.


### `$PRED` \ {#sec-block-pred}

__Syntax__: C++   
__Multiple allowed__: no   

Use `$PRED` to write a model without differential equations.  In this block, 
write all algebraic expressions for derived parameters, the response, and 
any other derived output quantities.  

For example:

```{c, eval=FALSE}
[ PARAM ] TVE0 = 100, AUC50 = 100, IMAX = 40, AUC = 0

[ PRED ]
double E0 = EVE0*exp(ETA(1));

double RESP = E0 - IMAX*AUC/(AUC50+AUC);
```

In this example, the entire model is written in the `$PRED` block.  It is an 
error to include the following blocks when `$PRED` is being used: `$MAIN`, 
`$TABLE`, `$PKMODEL`, `$ODE`, `$CMT`, `$INIT`.

See @sec-pred-data for additional information regarding data sets
in use with `$PRED` block.


### `$CAPTURE` \ {#sec-block-capture}

__Syntax__: text   
__Multiple allowed__: yes   
__Options__: `@annotated`, `@etas`

This is a block to identify variables that should be captured in the simulated
output. The `@etas` option is new with version 1.0.8.

For example:

```{c,eval=FALSE}
[ PARAM ] A = 1, B = 2

[ MAIN ]
double C = 3;
bool yes = true;

[ CAPTURE ] A B C yes
```

This construct will result in four additional columns in the simulated output 
with names `A`, `B`, `C`, and `yes`.  

Users can also rename captured variables by providing a 
`newname = oldname` specification. 

```{c, eval = FALSE}
$PARAM WT = 70, THETA1 = 2.2

$MAIN
double CL = THETA1*pow(WT/70,0.75)*exp(ETA(1));

$OMEGA 1

$CAPTURE WEIGHT = WT TVCL = THETA2 CL  ETA(1)
```

In this example, the names of the captured data items will be
`WEIGHT,TVCL,CL,ETA_1`.

Users can use the `capture` type to declare variables in `$MAIN` and `$TABLE`. 
`capture` types are really `double`s, but using that type will signal mrgsolve
to automatically capture that value.  For example:

```{c,eval=FALSE}
$PARAM VC = 300

$CMT CENT

$TABLE
capture DV = (CENT/VC);
```

Since we used type `capture` for `DV`, `DV` will show up as a column in the
simulated data.

Annotated example:

```{c,eval=FALSE}
$MAIN
double CLi = TVCL*exp(ECL);

$TABLE
double DV = (CENT/VC)*exp(PROP);

$CAPTURE @annotated
CLi : Individual clearance (L/hr)
DV  : Plasma concentration (mcg/ml)
```

Use the `@etas` option to specify an expression that evaluates to integers
which identify the `ETA` number(s) to capture. For example

```{c, eval = FALSE}
$OMEGA 0.1 0.2 0.3

$CAPTURE @etas 1:LAST
CL DV
```

Will capture all three `ETA`s into the simulation outputs. NONMEM-style
names will be used (e.g. `ETA1`). This is in contrast to the names that are 
generated when `ETA(1)` is captured; in this case, the parens are removed and 
replaced with `_` so the output name is `ETA_1`. mrgsolve will provide `LAST`
(and `last`) which will evaluate to the maximum number of `ETAs` in the problem
(in this case, 3). You can write any valid expression in `@etas`. When the 
expression is evaluated, it will be evaluated in the model environment created
through the `$ENV` block (@sec-block-env).


::: {.callout-tip}
In addition to listing model variables for output in the `$CAPTURE` block,
users can "dynamically" capture model variables through the `capture` 
argument to [`mread()`](https://mrgsolve.org/docs/reference/mread.html). These
may be model parameters (listed in `$PARAM`) or user-defined model variables.

To get a listing of available user-defined variables, you can use

```{r}
#| eval: false
mod <- house()
as.list(mod)$cpp_variables
```
:::


### `$OMEGA` \ {#sec-block-omega}

__Syntax__: text   
__Multiple allowed__: yes   
__Options__: `@annotated`, `@block`, `@correlation`, `@labels`, `@name`, `@object`, `@as_object`   

See `?modMATRIX` for more details about options for this block.

Use this block to enter variance/covariance matrices for subject-level random
effects drawn from multivariate normal distribution. All random effects are
assumed to have mean of 0.  Off diagonal elements for block matrices are assumed
to be correlation coefficients if the `@correlation` option is used (see below).

By default, a __diagonal__ matrix is assumed.  So:

```{c, eval=FALSE}
$OMEGA
1 2 3
```
will generate a 3x3 omega matrix.

A __block__ matrix may be entered by using `block = TRUE`. So:

```{c, eval = FALSE}
$OMEGA @block
0.1 0.02 0.3
```
will generate a 2x2 matrix with covariance 0.02.

A 2x2 matrix where the off-diagonal element is a correlation, not a covariance 
can be specified like this:

```{c,eval=FALSE}
$OMEGA @correlation
0.1 0.67 0.3
```

Here, the correlation is 0.67. mrgsolve will calculate the covariances and 
substitute these values.  The matrix will be stored and used with these 
covariances, not the correlation.

A name can be assigned to each matrix:
```{c, eval=FALSE}
$OMEGA @name PK @block
0.2 0.02 0.3

$OMEGA @name PD
0.1 0.2 0.3 0.5
```

to distinguish between multiple `$OMEGA` blocks and to facilitate updating 
later. The model in the preceding example will have two `$OMEGA` matrices: 
2x2 and 4x4.

Labels can be assigned which function as aliases for the different `ETAs`

```{c, eval = FALSE}
$OMEGA @block @labels ETA_CL ETA_V
0.1 0.05 0.2
```

The number of labels should match the number of rows (or columns) in the 
matrix.

Annotated example (diagonal matrix):

```{c, eval=FALSE}
$OMEGA @annotated
ECL: 0.09 : ETA on clearance
EVC: 0.19 : ETA on volume
EKA: 0.45 : ETA on absorption rate constant
```

Annotated example (block matrix):

```{c, eval=FALSE}
$OMEGA @annotated @block
ECL: 0.09 : ETA on clearance
EVC: 0.001 0.19 : ETA on volume
EKA: 0.001 0.001 0.45 : ETA on absorption rate constant
```

The `@object` and `@as_object` options should name or return a square numeric 
matrix. If `rownames` are included in the matrix, then they will be used to 
form labels for the realized ETAs. For example, we can initialize a very 
large `$OMEGA` matrix with 

```{c, eval = FALSE}
$OMEGA @as_object
n <- 20
lbl <- paste0("ETA_", LETTERS[1:n])
matrix(0, 20, 20, dimnames = list(lbl, lbl))
```

Note: this only initializes the matrix; you will (likely) need to update it 
with meaningful values after the model is loaded (@sec-matrix-chapter).


### `$SIGMA` \ {#sec-block-sigma}

__Syntax__: text   
__Multiple allowed__: yes   
__Options__: `@annotated`, `@block`, `@correlation`, `@labels`, `@name`, `@object`, `@as_object`  

See `?modMATRIX` for more details about options for this block.

Use this block to enter variance/covariance matrices for within-subject random
effects drawn from multivariate normal distribution. All random effects are
assumed to have mean of 0.  Off diagonal elements for block matrices are assumed
to be correlation coefficients if the `@correlation` option is used (see below).

The `@object` and `@as_object` options should name or return a square numeric 
matrix. If `rownames` are included in the matrix, then they will be used to 
form labels for the realized EPS values. 

The `$SIGMA` block functions like the `$OMEGA` block.  See `$OMEGA` for details. 

### `$SET` {#sec-block-set}

__Syntax__: R   
__Multiple allowed__: no   

Use this code block to set different options for the simulation. Use a 
`name = value` format, where `value` is evaluated by the `R` interpreter.  

Most of the options that can be entered in `$SET` are passed to `update`. 

For example:

```{c,eval=FALSE}
[ SET ] end = 240, delta = 0.5, req = "RESP"
```

Here, we set the simulation `end` time to 240, set the time difference between 
two adjacent time points to 0.25 time units, and request only the `RESP` 
compartment in the simulated output.


### `$GLOBAL`

__Syntax__: C++   
__Multiple allowed__: no   


The `$GLOBAL` block is for writing `C++` code that is outside of `$MAIN`, 
`$ODE`, and `$TABLE`.

There are no artificial limit on what sort of `C++` code can go in `$GLOBAL`.  
However there are two more-common uses:

1. Write `#define` preprocessor statements
1. Define global variables, usually variables other than `double`, `bool`, 
   `int` (see @sec-new-variables)

__Preprocessor directives__
Preprocessor  `#define` directives are direct substitutions that the `C++` 
preprocessor makes prior to compiling your code.  

For example:
```{c,eval=FALSE}
[ GLOBAL ]
#define CP (CENT/VC)
```
When this preprocessor directive is included, everywhere the preprocessor finds 
a `CP` token it will substitute `(CENT/VC)`.  Both `CENT` and `VC` must be 
defined and the ratio of `CENT` to `VC` will be calculated depending on whatever 
the current values are.  Notice that we included parentheses around `(CENT/VC)`.  
This makes sure the ratio between the two is taken first, before any other 
operations involving `CP`.  

__Declaring global variables__
Sometimes, you may wish to use global variables and have more control over how 
they get declared.  

```{c,eval=FALSE}
$GLOBAL
bool cure = false;
```
With this construct, the boolean variable `cure` is declared and defined right 
as the model is compiled.  

__Declare in bulk__

If you have a large number of variables to declare, you can do that in 
bulk in `$GLOBAL`. For example, we know we will need a long list of `double` 
precision variable for covariate modeling, we can declare them all at once like 
this:

```{c, eval = FALSE}
[ global ] 
double TVCL, TVV2, TVQ = 0, TVV3 = 0;

[ main ] 

TVCL = THETA1 * pow(WT / 70.0, 0.75);
TVV2 = THETA2 * WT / 70.0;
TVQ  = THETA3 * pow(WT / 70.0, 0.75);
TVV3 = THETA4 * WT / 70.0;
```

This isn't a terribly long list, but let's pretend it is to illustrate how to 
do this. You can also declare `int`, `bool` and other types like this. I've 
initialized the last two (`TVVQ` and `TVV3`) to illustrate how to do this. It's
good practice to do that but not necessary as long as everything gets
initialized to **something** before they are used.

### `$PKMODEL`

__Syntax__: R   
__Multiple allowed__: no   

This code block implements a one- or two-compartment PK model where the system 
is calculated by algebraic equations, not ODEs.  mrgsolve handles the 
calculations and an error is generated if both `$PKMODEL` and `$ODE` blocks are 
included in the same model specification file.

This is an options-only block. The user must specify the number of compartments
(1 or 2) to use in the model as well as whether or not to include a depot dosing
compartment.  See `?PKMODEL` for more details about this block, including
specific requirements for symbols that must be defined in the model
specification file.

The `$CMT` or `$INIT` block must also be included with an appropriate number of
compartments.  Compartment names, however, may be determined by the user.

Example:
```{c, eval=FALSE}
[ CMT ] GUT CENT PERIPH

[ PKMODEL ] ncmt=2, depot=TRUE
```

As of version `0.8.2`, we can alternatively specify the compartments
right in the `$PKMODEL` block:

```{c,eval=FALSE}
$PKMODEL cmt="GUT CENT PERIPH", depot = TRUE
```
Specifying three compartments with `depot=TRUE` implies `ncmt=2`.  Notice that a
separate `$CMT` block is not appropriate when `cmt` is specified in `$PKMODEL`.


### `$PLUGIN` \ {#sec-block-plugin}

__Syntax__: text   
__Multiple allowed__: yes   

Plugins are a way to add extensions to your mrgsolve model. Plugins can either 
link your model to external libraries (like `boost` or `Rcpp`) or they can 
open up access to additional functionality provided by mrgsolve itself. 

Plugins are listed and discussed in more detail in @sec-plugins.

__Usage__

To invoke a plugin, list the plugin name in the code block. For requesting
`Rcpp` headers in your model, call

```{c, eval = FALSE}
$PLUGIN Rcpp
```

__Available plugins__

The following plugins make additional mrgsolve-specific functionality available

- `mrgx`: extra `C++` functions (see below)
- `tad`: track time after dose in your model 
- `N_CMT`: get the number of a compartment

The following plugins give you a different model specification experience

- `autodec`: mrgsolve will find assignments in your model and automatically 
  declare them as `double`s
- `nm-vars`: a NONMEM look and feel for compartmental models; use `F1`, 
  `A(1)` and `DADT(1)` rather than `F_GUT`, `GUT` and `dxdt_GUT`

The following plugin lets you customize how your model is compiled

- `CXX11`: compile your model with `C++11` standard; this adds
  the compiler flag `-std=c++11`

The following plugins will let you link to external libraries

- `Rcpp`: include `Rcpp` headers int your model
- `BH`: include `boost` headers in your model 
- `RcppArmadillo`: include `Armadillo` headers in your model


Note that `Rcpp`, `RcppArmadillo` and `BH` only allow you to link to those
headers.  To take advantage of that, you will need to know how to use `Rcpp`,
`boost` etc.  For the `BH` plugin, no headers are included for you; you must
include the proper headers you want to use in `$GLOBAL`.


__mrgx__
This is a general collection of functions that we made available. 

Functions provided by `mrgx`:

  -  `T get<T>(std::string <pkgname>, std::string <objectname>)`
      - This gets an object of any Rcpp-representable type (`T`) from any package
  -  `T get<T>(std::string <objectname)`
      - This gets an object of any Rcpp-representable type (`T`) from `.GlobalEnv`
  -  `T get<T>(std::string <objectname>, databox& self)`
      - This gets an object of any Rcpp-representable type (`T`) from `$ENV`
  -  `double rnorm(double mean, double sd, double min, double max)`
      - Simulate one variate from a normal distribution that is between `min` and `max`
  -  `double rlognorm(double mean, double sd, double min, double max)`
      - Same as `mrgx::rnorm`, but the simulated value is passed to `exp` after simulating
  -  `Rcpp::Function mt_fun()` 
      - Returns `mrgsolve::mt_fun`; this is usually used when declaring a `R` function in `$GLOBAL`
      - Example: `Rcpp::Function print = mrgx::mt_fun();`


__IMPORTANT__ All of these functions are in the `mrgx` namespace.  So, in order
to call these functions you must include `mrgx::` namespace identifier to the
front of the function name.  For example, don't use `rnorm(50,20,40,140)`; use
`mrgx::rnorm(50,20,40,140)`.


#### Some examples

__Get a numeric vector from `$ENV`__
```{c,eval=FALSE}
[ PLUGIN ] Rcpp mrgx

[ ENV ]
x <- c(1,2,3,4,5)

[ GLOBAL ]
Rcpp::NumericVector x;

[ PREAMBLE ]
x = mrgx::get<Rcpp::NumericVector>("x", self);

```


__Get the `print` function from `package:base`__
```{c,eval=FALSE}
$PLUGIN Rcpp mrgx

$GLOBAL
Rcpp::Function print = mrgx::mt_fun();

$PREAMBLE
print = mrgx::get<Rcpp::Function>("base", "print");

$MAIN
print(self.rown);
```

Note that we declare the `print` in `$GLOBAL` and use the `mt_fun()` place holder.

__Simulate truncated normal variables__
This simulates a weight that has mean 80, standard deviation 20 and is greater
than 40 and less than 140.

```{c,eval=FALSE}
$PLUGIN Rcpp mrgx

$MAIN
if(NEWIND <=1) {
  double WT = mrgx::rnorm(80,20,40,140);
}
```

See also: @sec-block-preamble.

### `$THETA`  \ {#sec-block-theta}

__Syntax__: text   
__Multiple allowed__: yes   
__Options__: `@annotated`, `@name`, `@object`, `@as_object`

Use this code block as an efficient way to add to the parameter list where names
are determined by a prefix and a number.  By default, the prefix is `THETA` and
the number sequentially numbers the input values.

For example:

```{c,eval=FALSE}
[ THETA ]
0.1 0.2 0.3
```
is equivalent to
```{c,eval=FALSE}
$PARAM THETA1 = 0.1, THETA2 = 0.2, THETA3 = 0.3
```

Annotated example:

```{c,eval=FALSE}
$THETA @annotated
0.1 : Typical value of clearance (L/hr)
0.2 : Typical value of volume (L)
0.3 : Typical value of ka (1/hr)
```

To change the prefix, use `@name` option

```{c,eval = FALSE}
$THETA @name theta
0.1 0.2 0.3
```

would be equivalent to
```{c,eval = FALSE}
[ PARAM ] theta1 = 0.1, theta2 = 0.2, theta3 = 0.3
```

The `@object` and `@as_object` options should name or return an unnamed vector
of parameter values (names are ignored).

See also: @sec-block-param.


### `$NMXML`  \ {#sec-block-nmxml}

__Syntax__: R   
__Multiple allowed__: yes   

The `$NMXML` block lets you read and incorporate results from a NONMEM run into 
your mrgsolve model.


#### Syntax {.unnumbered}

If your model run is 1000 and the `.xml` file is located in `1000/1000.xml`, 
use this syntax

```{r}
#| eval: false
$NMXML 
run = 1000
project = "../model/nonmem"
root = "cppfile"
```

If your model run is 1000 and the `.xml` file is located in some other 
location use the `path` argument to specify the _complete_ path to the `.xml`
file

```{r}
#| eval: false
$NMXML 
path = "../<path>/<to>/<output>/<files>/1000.xml"
root = "cppfile"
```

This latter syntax may be needed when using PsN and output files are organized
in a way different from the previous example.


**Details**

From the NONMEM run, `THETA` will be imported into your
parameter list (see @sec-block-param and @sec-component-param), `OMEGA` will 
be captured as an `$OMEGA` block (@sec-block-omega) and `SIGMA` will be 
captured as a `$SIGMA` block (@sec-block-sigma).  Users may optionally omit 
any one of these from being imported.

`$NMXML` contains a `project` argument and a `run` argument.  By default, the 
estimates are read from from the file `project/run/run.xml`.  That is, it is 
assumed that there is a directory named `run` that is inside the `project` 
directory where `$NMXML` will find `run.xml`.  Your NONMEM run directories may
not be organized in a way that is compatible with this default.  In that case, 
you will need to provide the `path` argument, which should be the path to the 
`run.xml` file, either as a full path or as a path relative to the current 
working directory.

Once the model object is obtained, the path to the `xml` file that formed 
the source for imported parameters can be retrieved by coercing the 
model object to list and looking for `nm_import`:

```{r, eval = FALSE}
mod <- modlib("1005", compile = FALSE)

as.list(mod)$nm_import
```

For help on the arguments / options for `$NMXML`, please see the `?nmxml` 
help topic in your `R` session after loading the mrgsolve package.

__An example__

There is a NONMEM run embedded in the mrgsolve package 
```{r}
path <- file.path(path.package("mrgsolve"),"nonmem")
list.files(path, recursive=TRUE)

```

We can create a mrgsolve control stream that will import `THETA`, `OMEGA` and 
`SIGMA` from that run using the `$NMXML` code block.


```{mread, stem = "nmxml-1", annot = TRUE}
$NMXML
run = 1005
project = path
root = "cppfile"

olabels = c("ECL", "EVC", "EKA")
slabels = c("PROP", "ADD")

$MAIN
double CL = THETA1*exp(ECL);
double V2 = THETA2*exp(EVC);
double KA = THETA3*exp(EKA);
double Q = THETA4;
double V3 = THETA5;

$PKMODEL ncmt=2, depot=TRUE

$CMT GUT CENT PERIPH

$TABLE
double CP = (CENT/V2)*(1+PROP) + ADD/5;

$CAPTURE CP

$SET delta=4, end=96

```



__NOTE__: in order to use this code, we need to install the `xml2` package.

```{r, message = FALSE}
mod <- mread("inline/nmxml-1.cpp", compile = FALSE)

mod

param(mod)

revar(mod)

```

__root argument__: please use the `root = "cppfile"` argument going forward. 

As of mrgsolve `0.11.0`, we added an argument called `root` to `$NMXML` that 
tells mrgsolve the location where it should read the `xml` file from. The 
default behavior is the `"working"` directory. When this is the case, mrgsolve
assumes that the the `xml` file can be found relative to the "working" 
directory. The only other value that `root` can take is `"cppfile"`. When 
`root = "cppfile"`, then mrgsolve will look for the `xml` file in a directory 
that is relative to **where the model source code file is located**. Please
take a look at @sec-block-nmext for more discussion and examples. 

Starting with version 1.0.8, when the NONMEM run is found using the `run` +
`project` arguments, you can pass `run = "@cppstem"` and mrgsolve will look for
the NONMEM run number that matches the stem of the current mrgsolve file. For
example, if you are coding a model to match NONMEM run `101.ctl`, then call your
mrgsolve model file `101.cpp` or `101.mod`; just match the file stem. Then write

```{c, eval = FALSE, filename = "101.cpp"}
$NMXML 
run = "@cppstem"
project = <nonmem-run-directory>
```

mrgsolve will check the stem of the mrgsolve model file and  look for the NONMEM
`.xml` file in

```{r, eval = FALSE}
<nonmem-run-directory>/101/101.xml
```

Please see the `?nmxml` help topic for more information on arguments that can 
be passed to `$NMXML`.

See also: @sec-block-nmext.


### `$NMEXT`  \ {#sec-block-nmext}

__Syntax__: R   
__Multiple allowed__: yes   

Like `$NMXML`, `$NMEXT` allows the import of `$THETA`, `$OMEGA`, and `$SIGMA`
from your NONMEM run into your mrgsolve model, but the estimates are read from
the `.ext` file output. `$NMEXT` is able to import the NONMEM estimates much 
faster than `$NMXML` when loading from sampling based methods (mainly 
`METHOD=BAYES`). 

#### Syntax {.unnumbered}

If your model run is 1000 and the `.ext` file is located in `1000/1000.ext`, 
use this syntax

```{r}
#| eval: false
$NMEXT 
run = 1000
project = "../model/nonmem"
root = "cppfile"
```

If your model run is 1000 and the `.ext` file is located in some other 
location use the `path` argument to specify the _complete_ path to the `.ext`
file

```{r}
#| eval: false
$NMEXT 
path = "../<path>/<to>/<output>/<files>/1000.ext"
root = "cppfile"
```

This latter syntax may be needed when using PsN and output files are organized
in a way different from the previous example.

#### Details {.unnumbered}

Once the model object is obtained, the path to the `ext` file that formed 
the source for imported parameters can be retrieved by coercing the 
model object to list and looking for `nm_import`. 

```{r, eval = FALSE}
mod <- modlib("1005", compile = FALSE)

as.list(mod)$nm_import
```

See the `$NMXML` topic in @sec-block-nmxml for an example of equivalent 
functionality.


__IMPORTANT__: while `$NMEXT` works very similarly to `$NMXML`, there is one 
key difference between the two: when using `$NMEXT`, you will always get one 
`$OMEGA` matrix and one `$SIGMA` matrix, regardless of the block structure used
in the NONMEM control stream. When using `$NMXML`, you get the same block 
structure in the mrgsolve model as the NONMEM model.  For example: in the 
NONMEM control stream, `$OMEGA` is one 2x2 matrix and one 3x3 matrix. Importing
those estimates with `$NMEXT` will give you one 5x5 matrix, while importing 
the estimates with `$NMXML` will give you a list of one 2x2 and one 3x3 matrix.

__root argument__: please use the `root = "cppfile"` argument going forward. 

As of mrgsolve `0.11.0`, we added an argument called `root` to `$NMEXT` that 
tells mrgsolve the location where it should read the `ext` file from. The 
default behavior is the `"working"` directory. When this is the case, mrgsolve
assumes that the the `ext` file can be found relative to the "working" 
directory. The only other value that `root` can take is `"cppfile"`. When 
`root = "cppfile"`, then mrgsolve will look for the `ext` file in a directory 
that is relative to **where the model source code file is located**. 

We recommend that users start using the `root` argument and set it to "cppfile". 
This will eventually become the default. For example:

```{c, eval = FALSE}
$NMEXT
run = 1005
project = "../../model/pk"
root = "cppfile"
```

This tells mrgsolve to find the nonmem run back two directories (`../../`) 
and then into `model --> pk --> 1005` relative to where the mrgsolve model 
file is located. This is in contrast to the previous expected behavior that 
the path should be relative to the current working directory.

We also note here that when users pass an absolute path, the relative path 
doesn't matter at all. Users are free to use pathing tools to generate the 
absolute path to the project.  For example, the `here::here()` function can 
be used like this

```{c, eval = FALSE}
$NMEXT
run = 1005
project = here::here("model/pk")
```

When used in the proper context, `here()` will generate an absolute path from 
your projects root directory to the nonmem project directory. Please refer to 
the `here::here()` documentation for proper use of this function. 

Starting with version 1.0.8, when the NONMEM run is found using the `run` +
`project` arguments, you can pass `run = "@cppstem"` and mrgsolve will look for
the NONMEM run number that matches the stem of the current mrgsolve file. 
See details and example in @sec-block-nmxml.

See the `?nmext` R help topic for arguments that can be passed to `$NMEXT` block.
Notably, the user can select the function to read in the ext file.By default, 
mrgsolve will try to load data.table and use the `fread` function. If 
data.table can't be loaded, then mrgsolve will use `utils::read.table`.

See also: @sec-block-nmxml.

### `$INCLUDE`  \ {#sec-block-include}

__Syntax__: text   
__Multiple allowed__: no   

To include your own header file(s) in a model use `$INCLUDE`
```{c,eval=FALSE}
$INCLUDE
mystuff.h
otherstuff.h
```

or

```{c,eval=FALSE}
$INCLUDE
mystuff.h,  otherstuff.h
```

mrgsolve will insert proper `#include` preprocessor directives into the `C++` 
code that gets compiled.

__Requirements__

  - All header files listed in `$INCLUDE` are assumed (expected) to be under 
    the `project` directory; don't use `$INCLUDE` for header files that are in 
    any other location
  - An error is generated if the header file does not exist
  - An error is generated if any quotation marks are found in the file name 
    (don't use quotes around the file name; mrgsolve will take care of that)
  - A warning is issued if the header file does not end in `.h`
  - When the header file is changed (MD5 checksum changes), the model will be 
    forced to be rebuilt (recompiled) when `mread()` or `mcode()` (but not 
    `mread_cache()` or `mcode_cache()`) is called; this feature is only available 
    for header files listed in `$INCLUDE` (see below)
  - Do not use `$INCLUDE` to include `Rcpp`, `boost`, `RcppArmadillo` or 
    `RcppEigen` headers; use the appropriate `$PLUGIN` instead

For applications that don't fit into the requirements listed above, users can 
always include header files in the model in `$GLOBAL` like this:

```{c,eval=FALSE}
$GLOBAL
#include "/Users/me/libs/mystuff.h"
```

But be careful when doing this: if there are changes to `mystuff.h` but not to 
any other part of the model specification, the model may not be fully compiled 
when calling `mread`.  In this case, always use `preclean=TRUE` argument to 
`mread` to force the model to be built when calling `mread`.

### `$ENV` \ {#sec-block-env}

__Syntax__: R   
__Multiple allowed__: no   

This block is all `R` code (just as you would code in a stand-alone `R` script. 
The code is parsed and evaluated into a new environment when the model is 
compiled.  

For example:

```{c,eval=FALSE}
$ENV

Sigma <- cmat(1,0.6,2)

mu <- c(2,4)

cama <- function(mod) {
  mod %>%
    ev(amt=100, ii=12, addl=10) %>% 
    mrgsim(obsonly=TRUE,end=120)
}

```

Objects inside `$ENV` can be utilized in different C++ functions 
(see @sec-block-plugin) or other parts of the simulation process.

Objects inside `$ENV` will also be made available when parsing different
parts of the model specification file. For example, we can initialize an
`$OMEGA` block with the `Sigma`  object in the `$ENV` block example above:

```{c, eval = FALSE}
$OMEGA @object Sigma
```


Another example is providing constants or functions to simplify code in the `$PARAM` block

```{c, eval = FALSE}
$ENV
rescale <- 1/1000 

convert <- function(x) x * 5 / 166.2

$PARAM
p = 500 * rescale

q  = convert(2.51)
```

You can also make your model specification file depend on an R package; 
for example, if you want to use the `here` package to render a path to 
your NONMEM model

```{c, eval = FALSE}
$ENV
if(!requireNamespace("here")) {
  stop("the `here` package is required to load this model.")  
}

$NMXML
run = 101
project = here::here("model/nonmem")
```

You can access the model environment using the `env_get()` function: 

```{c, eval = FALSE}
env_get(mod)
```

or directly in the `envir` slot

```{c, eval = FALSE}
mod@envir
```

See also `env_ls()` to list the names of variables in the model environment, 
`env_update()` to change the values of objects in the model environment, 
or `env_eval` to re-evaluate the model environment code.   


## Variables and Macros {#sec-variables-macros}

This section describes some macros and internal variables that can be used in
model specification files.  In the following section, we adopt the convention
that `CMT` stands for a compartment in the model.

__IMPORTANT NOTE__:

It should be clear from the usage examples which variables can be set by the
user and which are to be read or checked.  **All internal variables are
pre-defined and pre-initialized by mrgsolve. _The user should never try to
declare an internal variable; this will always result in an compile-time
error_**.


### `ID`
The current subject identifier. `ID` is an alias for `self.id`.

### `TIME`
Gives the time in the current data set record.  This is usually only used in 
`$MAIN` or `$TABLE`.  `TIME` is an alias for `self.time`.  Contrast with 
`SOLVERTIME`.

### `SOLVERTIME`
Gives the time of the current timestep taken by the solver.  *This is can only 
be used in `$ODE`*.  Contrast with `TIME`.

Starting with mrgsolve version 1.0.0, the variable `T` is provided as a synonym 
to `SOLVERTIME` inside `$ODE`  when you invoke the `nm-vars` plugin
(@sec-plugin-nm-vars).

### `T` 
Only provided when the `nm-vars` plugin is invoked; see `SOLVERTIME`.

### `EVID`
`EVID` is an event id indicator.  mrgsolve recognized the following event IDs:

* 0 = an observation record
* 1 = a bolus or infusion dose
* 2 = other type event, with solver reset
* 3 = system reset
* 4 = system reset and dose
* 8 = replace

`EVID` is an alias for `self.evid`.

### `CMT`

`CMT` is the current compartment number.  In this case, `CMT` is used literally
(not a stand-in for the name of a compartment).   For example:

```{c, eval=FALSE}
$CMT GUT CENT PERIPH

$MAIN

if(CMT==2) {
  // ....
}
```

In the example, `GUT`, `CENT` and `PERIPH` are the amounts in the respective
compartments and `CMT` refers to the value of `CMT` in the data record / 
data set.

`CMT` is an alias for `self.cmt`.

### `AMT`

`AMT` is the current value of dose amount.

`AMT` is an alias for `self.amt`.


### `NEWIND`
`NEWIND` is a new individual indicator, taking the following values:

* 0 for the first event record of the data set
* 1 for the first event record of a subsequent individual
* 2 for subsequent event record for an individual

For example:
```{c, eval=FALSE}
[ GLOBAL ]
int counter = 0;

[ MAIN ]
if(NEWIND <=1) {
  counter  = 0;
}
```

`NEWIND` is an alias for `self.newind`.

### `SS_ADVANCE`

This is a `bool` data item (either `true` or `false`) which is always `false` 
unless mrgsolve is currently advancing the system to steady state
(@sec-steady-state); then it will be `true`. This variable is only 
available in the `$ODE` (or `$DES`) block. 

Use this variable to modify the calculation of your model differential equations
while the system is advancing to steady state. One use case is when you have
an accumulation compartment to calculate `AUC` but you want to halt accumulation
while the system is working toward steady state

```{c, eval = FALSE}
$ODE
dxdt_CENT = -(CL/V) * CENT;
dxdt_AUC = CENT/V;
if(SS_ADVANCE) dxdt_AUC = 0;
```


### `simeta()`

The `simeta()` function can be used to re-simulate ETA values. For example, 

```{c, eval = FALSE}
$MAIN
simeta();
```

will re-simulate all `ETA(n)`. 

Starting with mrgsolve version 1.6.1, `simeta()` is available in `$MAIN`
(`$PK)`, `$TABLE` (`$ERROR`), `$EVENT` and `$PRED` blocks. 


### `simeps()`

`simeps()` works like `simeta()`, but all `EPS(n)` values are re-simulated rather
than the ETA values.


### `self` \ {#sec-self-object}

`self` is an object that gets passed to your model function that contains 
data items and functions that can be called. It is a `struct`  (see
the source code 
[here](https://github.com/metrumresearchgroup/mrgsolve/blob/master/inst/base/mrgsolv.h)). 
A partial list of members are documented in the following sections. 

`self` functions include

- `self.tad()`
- `self.mtime()`
- `self.mevent()`
- `self.stop()`
- `self.stop_id()`
- `self.stop_id_cf()`

`self` data members include

- `self.id` (`ID`)
- `self.amt` (`AMT`)
- `self.cmt` (`CMT`)
- `self.evid` (`EVID`)
- `self.time` (`TIME`)
- `self.newind` (`NEWIND`)
- `self.nid` 
- `self.idn`
- `self.rown`
- `self.nrow`
- `self.envir` (see note below)

This information is provided for transparency and is not exhaustive. 
We provide an interface to these data items through pre-processor directives
with simpler names (e.g. `EVID` will translate to `self.evid`) and users are 
encouraged to use the simpler, API name. 

### `self.time` {#sec-self.time}

The current data set `TIME`.

### `self.cmt`
The current compartment number regardless of whether it was given as `cmt` or 
`CMT` in the data set.  There is no alias for `self.cmt`.

For example:

```{c,eval=FALSE}
$TABLE

double DV = CENT/VC + EPS(1);

if(self.cmt==3) DV = RESPOSE + EPS(2);

```

### `self.amt`
The current `amt` value regardless of whether it was given as `amt` or `AMT` 
in the data set.  There is no alias for `self.amt`. 

```{c, eval=FALSE}
[ PREAMBLE ]
double last_dose = 0;

[ MAIN ]

if(EVID==1) {
  last_dose = self.amt;  
}

```

### `self.nid`
The number of IDs in the data set.

### `self.idn`
The current ID number.  Numbers start at 0 and increase by one to `self.nid-1`. 
So if you want to test for the last ID in the output data set, you would write:

```{c, eval = FALSE}
[ table ]

if(self.idn == (self.nid-1)) {
  // do something ....
}
```

### `self.nrow`
The total number of rows in the output data set.

### `self.rown`
The current row number.  Numbers start at 0 and increase by one to `self.rown-1`. 
So if you want to test for the last row in the output data set, you would write:

```{c, eval = FALSE}
[ table ]

if(self.rown == (self.nrow-1)) {
  // do something ....
}
```

### `self.envir`

This item is a null pointer that can be cast to an `Rcpp::Environment` only 
when the `Rcpp` plugin is invoked (see @sec-plugin-rcpp).

For example

```{c, eval = FALSE}
[ plugin ] Rcpp mrgx
  
[ preamble ] 
Rcpp::Environment env = mrgx::get_envir(self);
```



### `self.tad()` \ {#sec-self.tad}

This is a function that calculates and returns the time after the most recent
dose event record (any record with EVID equal to 1 or 4). `self.tad()` will 
return `-1` when it is called before the first dose record within an individual
(`NEWIND <= 1`). 

This function should be called in `$MAIN` every time `$MAIN` is called.

For example:

```{c, eval = FALSE}
$MAIN

double TAD = self.tad();
```

Do not make this calculation depend on any test or condition; it must be called
every time `$MAIN` is called in order to see every dose.


### `self.mtime(<time>)` \ {#sec-self.mtime}

This is a function that creates a modeled even time.  The only argument is a 
numeric time in the future for the current individual that indicates when a 
discontinuity should be included in the simulation.  

When when `self.mtime()` is called, a new record is added to the current 
individual's record set with `EVID` set to 2.  This means that when the system
advances to this record (this time), the differential equation solver will 
reset and restart, creating the discontinuity.  The function returns the time of
this event so the user can work with it in subsequent code. 

For example, 

```{c, eval = FALSE}
$PARAM change_point = 5.13;

$MAIN

double KA = 1.1;

double mt1 = self.mtime(change_point);

if(TIME >= mt1) KA = 2.1;  

```

### `self.mevent(<time>, <evid>)`

Related to `self.mtime()` (@sec-self.mtime), except you can set a 
specific `EVID` for this intervention and you track the change time via `EVID`
rather than the time. You'd use this if you want to anchor 
several events in the future and be sure you can distinguish between them. 
For example


```{c, eval = FALSE}
[ main ] 
self.mevent(change_point1, 33);
self.mevent(change_point2, 34);
```

Now, you can look for `EVID==33` to know you have hit `change_point1` or 
`EVID==34` to indicate that you have hit `change_point2`. Notice that 
`self.mevent()` doesn't return any value. You will know that you've hit the 
change time by checking `EVID`.

### `self.stop()`

This `self` function is available to be called from `$PREAMBLE`, `$MAIN`, 
and `$TABLE`. When this function is called, the entire problem is stopped upon
processing the next simulation record. 

This might be called when something really bad happened and you just want to 
stop the simulation with an error.

### `self.stop_id()`

This `self` function is available to be called from `$PREAMBLE`, `$MAIN`, 
and `$TABLE`. When this function is called, processing of the current individual
is stopped and missing values (`NA_real`) are filled in for remaining compartment'
and capture outputs. 

This might be called when some condition is reached in the current individual
that indicates either that the rest of the outputs are inconsequential or there
was a problem with this particular individual.

See also `self.stop_id()` and `self.stop()`.

### `self.stop_id_cf()`

This `self` function is available to be called from `$PREAMBLE`, `$MAIN`, 
and `$TABLE`. When this function is called, processing of the current individual
is stopped and current values are carried forward (`cf`) for the remaining
output records for that individual. 

This might be called when some condition is reached in the current individual
that indicates either that the rest of the outputs are inconsequential or there
was a problem with this particular individual.

See also `self.stop_id_cf()` and `self.stop()`.


### `THETA(n)`

`THETA(n)` is ubiquitous in NONMEM control streams, representing estimated 
fixed effect parameters. mrgsolve does not attach any special meaning to 
`THETA(n)` but starting with version 1.0.0 it will translate `THETA(n)` to 
`THETAn`.  This is useful when importing NONMEM estimates using `$NMXML`
or `$NMEXT`, where `THETA`s are imported as parameters with names 
`THETA1`, `THETA2`, `THETA3` and so on (see @sec-block-nmxml and
@sec-block-nmext). `THETA(n)` is simply a NONMEM-oriented style for referring
to `THETAn`.

### `ETA(n)`
`ETA(n)` is the value of the subject-level variate drawn from the model 
`OMEGA` matrix.  `ETA(1)` through `ETA(25)` have default values of zero so 
they may be used in a model even if appropriate `OMEGA` matrices have not been
provided. 

For example:
```{c, eval=FALSE}
$OMEGA
1 2 3

$MAIN
double CL = TVCL*exp(ETA(1));
double VC = TVVC*exp(ETA(2));
double KA = TVKA*exp(ETA(3));
```
Here, we have a 3x3 `OMEGA` matrix. `ETA(1)`, `ETA(2)`, and `ETA(3)` will be
populated with variates drawn from this matrix.  `ETA(4)` through `ETA(25)` 
will be populated with zero. 

### `EPS(n)`
`EPS(n)` holds the current value of the observation-level random variates drawn
from `SIGMA`.  The basic setup is the same as detailed in `ETA(n)`.

Example:

```{c,eval=FALSE}
[ CMT ] CENT

[ PARAM ] CL=1, VC=20

[ SIGMA ] 
25 0.0025

[ TABLE ]
double DV = (CENT/VC)*(1+EPS(2)) + EPS(1);

```

### `SIGMA(n)`

Starting with version 1.0.8, users have read-only access to on-diagonal elements
of $\Sigma$ through the `SIGMA()` macro. For example

```{c, eval = FALSE}
$SIGMA 0.1 12

$TABLE
double STD = sqrt(SIGMA(1) + pow(F,2)*SIGMA(2));
```


### `table(<name>) `
This macro has been deprecated.  Users should __not__ use code like this:

```{c,eval=FALSE}
[ TABLE ]
table(CP) = CENT/VC;
```

But rather this:
```{c,eval=FALSE}
$TABLE 
double CP = CENT/VC;

$CAPTURE CP
```

See: @sec-block-table and also @sec-block-capture.

### `END_OF_INFUSION` \ {#sec-spec-end-of-infusion}

This is a boolean variable indicating if the current record
terminates an infusion. 

```{c}
/*| eval: false */
$ERROR

if(END_OF_INFUSION) {
  // Do something ... 
}
```

### `F_CMT` \ {#sec-spec-f}
For the `CMT` compartment, sets the bioavailability fraction for that 
compartment. 

Example:
```{c,eval=FALSE}
$MAIN
F_CENT = 0.7;
```

### `Fn`
Only available when the `nm-vars` plugin (@sec-plugin-nm-vars) is invoked; 
see `F_CMT`.


### `R_CMT` \ {#sec-spec-r}
For the `CMT` compartment, sets the infusion rate for that compartment.  The
infusion rate is only set via `R_CMT` when `rate` in the data set or event 
object is set to `-1`. 

Example:
```{c,eval=FALSE}
$MAIN
R_CENT = 100;
```

### `Rn`
Only available when the `nm-vars` plugin (@sec-plugin-nm-vars) is invoked; 
see `R_CMT`.


### `D_CMT` \ {#sec-spec-d}
For the `CMT` compartment, sets the infusion duration for that compartment.  
The infusion duration is only set via `D_CMT` when `rate` in the data set or 
event object is set to `-2`. 

Example:
```{c,eval=FALSE}
$MAIN
D_CENT = 2;
```

Starting with mrgsolve version 1.6.1, `D_CMT` is available in `$MAIN`
(`$PK)`, `$TABLE` (`$ERROR`), `$EVENT` and `$PRED` blocks. 

### `Dn`
Only available when the `nm-vars` plugin (@sec-plugin-nm-vars) is invoked; 
see `D_CMT`.

### `ALAG_CMT` \ {#sec-spec-alag}
For the `CMT` compartment, sets the lag time for doses into that compartment. 

Example:
```{c,eval=FALSE}
$MAIN
ALAG_GUT = 0.25;
```

Starting with mrgsolve version 1.6.1, `ALAG_CMT` is available in `$MAIN`
(`$PK)`, `$TABLE` (`$ERROR`), `$EVENT` and `$PRED` blocks. 


### `ALAGn`
Only available when the `nm-vars` plugin (@sec-plugin-nm-vars) is invoked; 
see `ALAG_CMT`.


### NONMEM-like syntax

There is a NONMEM-like syntax that lets you write `F1` rather than `F_GUT`, 
`D2` rather than `D_CENT` as well as other variables that are commonly used
in a NONMEM control stream. To make this syntax available, the `nm-vars` 
plugin must be invoked. 

See @sec-plugin-nm-vars on the `nm-vars` plugin for details. 

## Reserved words

Reserved words cannot be used as names for parameters, compartments or other 
derived variables in the model. **Note** that some of these words are "reserved"
for you to use in your data set. 

```{r, comment = " ", eval = FALSE, code = mrgsolve:::Reserved}

```


Other reserved words depend on the compartment names in your model. For 
example, if you have a compartment called `CENT` in the model, then the 
the following will be reserved


```{r, eval = FALSE}
F_CENT
R_CENT
D_CENT 
ALAG_CENT
N_CENT
```


## Derive new variables {#sec-new-variables}

New `C++` variables may be derived in `$GLOBAL`, `$PREAMBLE` `$MAIN`, `$ODE` and
`$TABLE`. Because these are `C++` variables, the type of variable being used
must be declared. For the vast majority of applications, the `double` type is
used (double-precision numeric value).

```{c,eval=FALSE}
$MAIN

double CL = TVCL*exp(ETA(1));
```

We want `CL` to be a numeric value, so we use `double`.  To initialize a `boolean` 
variable (true / false), write
```{c,eval=FALSE}
$MAIN
bool cure = false;
```

### Special handling for `double`, `int`, `bool`

When variables of the type `double`, `int`, and `bool` are declared and
initialized in `$PREAMBLE`, `$MAIN`, `$ODE`, `$TABLE` or `$EVENT`, mrgsolve will 
detect those declarations, and modify the code so that the variables are actually
declared once in `$GLOBAL` not in `$MAIN`, `$ODE`, or `$TABLE`.  This is done so
that variables declared in one code block (e.g. `$MAIN`) can be read and
modified in another code block (e.g. `$TABLE`).

For example, in the following code:
```{c,eval=FALSE}
$MAIN
double CL = TVCL*exp(ETA(1));
```

a double-precision numeric variable is created (`CL`) in the `$MAIN` block. 
When mrgsolve parses the model file, this code gets translated to

```{c,eval=FALSE}
$GLOBAL
namespace {
  double CL;
}

$MAIN
CL = TVCL*exp(ETA(1));
```

That is, `CL` is declared in `$GLOBAL` in an unnamed namespace so that variables
like this are global variables within the model file only.  

This way, we can still read the `CL` variable in `$TABLE`:

```{c,eval=FALSE}
$MAIN
double CL = TVCL*exp(ETA(1));
double V2 = TVV2*exp(ETA(2));

$TABLE
double KE = CL/V2;

$CAPTURE KE
```

To declare a variable that is local to a particular code block:

```{c,eval=FALSE}
$MAIN

local_double CL = TVCL*exp(ETA(1));
```

The `local_double` type is still just a double-precision variable.  The
difference is that it is protected from this re-declaration process and the
variable will be local to (in this case) the `$MAIN` block.

Similarly, use `local_int` or `local_bool` to declare a variable that will _not_
be made global by mrgsolve. One common application of `local_int` is when 
you need to write a for-loop

```{c,eval=FALSE}
$TABLE 
for(local_int i = 0; i < 9; ++i) {}. 
  // do something with i
}
```

Here, we we want `i` declared locally for use in that for-loop and we use 
`local_int` to make that happen.

### Initializing and resetting 

mrgsolve takes a hands-off approach to variables created by the user; it will
not initialize or reset variables from record to record or from individual to
individual. Failing to properly initialize or reset user variables can have 
unintended consequences, especially when conditional logic is used. 

Generally, a default value for all variables should be set prior to 
conditional assignment based on some test condition. 


### Using other types globally

As we noted in the previous section, `double`, `int`, and `bool` are processed
in a special way so that they are by default global to the file. Many times 
we want to work with other variable types in a global manner.  Whenever you 
want a data structure to be accessible across functions (e.g. `$MAIN`, `$TABLE`, 
etc.) they should be declared in `$GLOBAL`, optionally in an unnamed namespace.

For example:

```{c,eval=FALSE}
[ GLOBAL ]
std::vector<double> myvec;
```

or

```{c,eval=FALSE}
[ GLOBAL ]
namespace {
  std::vector<double> myvec;
}
```

In case that object needs some configuration prior to starting the problem, 
use `$PREAMBLE` to do that work

```{c, eval=FALSE}
[ GLOBAL ]
std::vector<double> myvec;

[ PREAMBLE ] 
myvec.assign(3,1);
```

### Automatic declaration

There is a plugin you can invoke which will instruct mrgsolve to scan your
code and discover variables to be declared with type `double`. See the 
autodec plugin in @sec-plugin-autodec for details.

## Random number generation

Users can simulate random numbers inside a model file using functions that 
are similar to the functions you'd normally use in R (e.g. `rnorm()` and 
`runif()`). This functionality is provided by `Rcpp` and therefore requires
using the `Rcpp` plugin (see @sec-plugins and @sec-plugin-rcpp). 

`Rcpp` provides these functions inside the `R` namespace so you will have to 
prefix the function call with `R::`. 

As an example, to make a draw from Uniform (0,1)

```{c, eval = FALSE}
[ plugin ] Rcpp
  
[ error ] 
double draw = R::runif(0,1);
```

Note that the `0` gets used as `min` and the `1` gets used as `max`;  
we didn't pass `n` here and  `draw` is a single number (not a vector like you'd 
get from (`runif(100, 0, 1)` on your R console). So in general, these functions 
work like their `R` counterparts, but without the `n` argument.

Another example showing how to draw from binomial distribution with probability
0.5

```{c, eval = FALSE}
[ plugin ] Rcpp
  
[ error ] 
double draw = R::rbinom(1, 0.5);
```

Here, the `1` is used as `size` (not `n`) and `0.5` is used as `prob`. 


Other helpful functions could be `R::rnorm()` or `R::rlnorm()` but you can 
call any of the `r` functions as well as the corresponding `dpq` functions
through this `R` namespace. 

Documentation, including functions to call and arguments can be found in the 
Rcpp API docs

http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html


## Examples

The following sections show example model specification.  The intention is to 
show how the different blocks, macros and variables can work together to make 
a functional model.  Some models are given purely for illustrative purpose and 
may not be particularly useful in application.

### Simple PK model

Notes:

* Basic PK parameters are declared in `$PARAM`; every parameter needs to be
  assigned a value
* Two compartments `GUT` and `CENT` are declared in `$CMT`; using `$CMT` 
  assumes that both compartments start with 0 mass
* Because we declared `GUT` and `CENT` as compartments, we write `dxdt_` 
  equations for both in `$ODE`
* In `$ODE`, we refer to parameters (`CL/VC/KA`) and the amounts in each
  compartment at any particular time (`GUT` and `CENT`)
* `$ODE` should be `C++` code; each line ends in `;`
* We derive a variable called `CP` in `$TABLE` that has type `capture`; mrgsolve
  will enter the `CP` name into the `$CAPTURE` block list

```{c, eval=FALSE}
$PARAM CL = 1, VC = 30, KA = 1.3

$CMT GUT CENT

$ODE

dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - (CL/VC)*CENT;

$TABLE
capture CP = CENT/VC;

```

This model can also be written without differential equations

```{c, eval=FALSE}
[ PARAM ] CL = 1, VC = 30, KA = 1.3

[ PKMODEL ] cmt = "CMT GUT CENT", depot = TRUE

$TABLE
capture CP = CENT/VC;

```
### PK/PD model

Notes:

* We use a preprocessor `#define` directive in `$GLOBAL`; everywhere in the 
  model where a `CP` token is found, the expression `(CENT/VC)` ... with 
parentheses ... is inserted
* We write the initial value for the `RESP` compartment in `$MAIN` as a 
  function of two parameters `KIN/KOUT`
* A new variable - `INH`- is declared and used in `$ODE`
* Since `CP` is defined as `CENT/VC`, we can "capture" that name/value in
`$CAPTURE`
* Both `$MAIN` and `$ODE` are `C++` code blocks; don't forget to add the `;` 
at the end of each statement

```{c, eval=FALSE}
$PARAM CL = 1, VC = 30, KA = 1.3
KIN = 100, KOUT = 2, IC50 = 2

$GLOBAL
#define CP (CENT/VC)

$CMT GUT CENT RESP

$MAIN
RESP_0 = KIN/KOUT;

$ODE

double INH = CP/(IC50+CP);

dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - (CL/VC)*CENT;
dxdt_RESP = KIN*(1-INH) - KOUT*RESP;

$CAPTURE CP
```

### Population PK model with covariates and IOV

Notes:

* Use `$SET` to set the simulation time grid from 0 to 240 by 0.1
* There are two `$OMEGA` matrices; we name them `IIV` and `IOV`
* The IIV "etas" are labeled as `ECL/EVC/EKA`; these are aliases to
  `ETA(1)/ETA(2)/ETA(3)`.  The `IOV` matrix is unlabeled; we must refer to
  `ETA(4)/ETA(5)` for this
* Because `ETA(1)` and `ETA(2)` are labeled, we can "capture" them as `ECL` and
  `EVC`
* We added zeros for both `$OMEGA` matrices; all the etas will be zero until we
  populate those matrices (@sec-matrix-update)


```{c,eval=FALSE}
$PARAM TVCL = 1.3, TVVC=28, TVKA=0.6, WT=70, OCC=1

$SET delta=0.1, end=240

$CMT GUT CENT

$MAIN

double IOV = IOV1
if(OCC==2) IOV = IOV2;

double CLi = exp(log(TVCL) + 0.75*log(WT/70) + ECL + IOV);
double VCi = exp(log(TVVC) + EVC);
double KAi = exp(log(TVKA) + EKA);

$OMEGA @name IIV @labels ECL EVC EKA
0 0 0
$OMEGA @name IOV @labels IOV1 IOV2
0 0

$SIGMA 0

$ODE
dxdt_GUT = -KAi*GUT;
dxdt_CENT = KAi*GUT - (CLi/VCi)*CENT;

$TABLE
capture CP = CENT/VCi;

$CAPTURE IOV ECL EVC
```