93  Stable Routines for Estimation Methods and Automated Protection Against Floating Point Exceptions (nm74)

Introduction to NONMEM 7.6.0

As of nm74, a series of routines are available that protect against domain violations, divide by zero, and floating point overflows, so you won’t need to insert the extra protective code lines described in the previous section. Each of these routines start with the letter P, followed by the name of the mathematical operation they are to perform. For example, PLOG is the protective code routine that performs the LOG operation. In addition, there are first derivative (such as PLOGD1), and second derivative (such as PLOGD2) companion routines available which NONMEM uses for analytical derivatives. The source code of these routines are available in ..\source\PROTECT.f90. If you wish to modify their behavior, then copy PROTECT.f90 to your run directory, rename and modify it, such as PROTECTB.f90, then refer to this modified code with

$SUBROUTINES OTHER=PROTECTB.f90

Make sure you modify the companion *D1 and *D2 routines appropriately as well.

The following protective code routines are available:

For all routines, if X=not a number, X is converted to machine precision value, which is about 1.0E-15, before performing an operation on it. If X>INFNTY (where INFNTY is approximately 1.0E+154), then X is converted to INFNTY before an operation is performed on it.

PLOG(x): returns LOG of x. If x<TINYZ, where TINYZ is approximately 2.2E-308, then -SQINFN=-1.1579E+077 is returned.

PLOG10(x): returns LOG10 of x. If x<TINYZ, then -SQINFN/LOG(10) is returned.

PSQRT(x): returns SQRT of x. If x<=0.0d+00, then 0 is returned.

PEXP(x): returns EXP of x. If x>350.0, then PEXP(350.0) is returned (avoids overflow).

PDZ(x): returns 1/x. Protects against divide by zero. If abs(x)<SMALLZ, then 1/SMALLZ is returned.

PZR(x): returns x. protects against zero. If abs(x)<SMALLZ, then SMALLZ is returned.

PNP(x): returns x. Protects against non-positive. If X<SMALLZ, then SMALLZ is returned.

PHE(x): returns x. Protects against high exponent. If X>350, then 350 is returned. Thus PEXP(x)=EXP(PHE(x)).

PNG(x): returns x. Protects against negative. If X<0.0, then 0.0 is returned.

PTAN(x): returns tan(x). Protects against returning infinity on inputs near pi/2.

PATAN(x): returns atan(x). Protects against large intputs.

PACOS, PASIN, returns acos(x), asin(x), respectively. If |X| is between 1.0 and 1+10**(-08), then x is submitted as 1 or -1. So, “dirty ones” are cleaned up, but values clearly beyond 1 are allowed to trip up the function, so the user is aware of the logical error in the code, and fix the issue.

Instead of replacing various operations with protected code operations by hand, you can ask NMTRAN to automatically convert your code to protected code with the following statement:

$ABBR PROTECT

NMTRAN will automatically replace all LOG (or DLOG) with PLOG, EXP (or DEXP) with PEXP, SQRT (or DSQRT) with PSQRT, / operations with *PDZ(), and B**E operations with PEXP(E*PLOG(B)).

When you use $ABBR PROTECT, you will find a considerable improvement in estimation stability, regardless of estimation method used.