Compartmental Analysis

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: IMNLPEX3                                            */
/*   TITLE: Compartmental Analysis                              */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS:                                                     */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: saswmh/sasghg               UPDATE: frwick 5/2015   */
/*     REF: Technical Report Examples                           */
/*    MISC:                                                     */
/*                                                              */
/****************************************************************/

data tetra;
   input t c @@;
   datalines;
 1 0.7   2 1.2   3 1.4   4 1.4   6 1.1
 8 0.8  10 0.6  12 0.5  16 0.3
;

proc iml;
use tetra;
read all into tetra;
close tetra;

start f(theta) global(thmtrx,t,h,tetra,eps);
   thmtrx = ( -theta[1] || 0 )     //
            (  theta[1] || -theta[2] );
   c = theta[3]//0 ;
   t = 0 // tetra[,1];
   call ode( r1, "der",c , t, h) j="jac" eps=eps;
   f = ssq((r1[2,])`-tetra[,2]);
   return(f);
finish;

start der(t,x) global(thmtrx);
   y = thmtrx*x;
   return(y);
finish;

start jac(t,x) global(thmtrx);
   y = thmtrx;
   return(y);
finish;

h      = {1.e-14 1. 1.e-5};
opt    = {0 2 0 1 };
tc     = repeat(.,1,12);
tc[1]  = 100;
tc[7]  = 1.e-8;
par    = { 1.e-13 . 1.e-10 . . . . };
con    = j(1,3,0.);
itheta = { .1  .3  10};
eps    = 1.e-11;

call nlpdd(rc,rx,"f",itheta) blc=con opt=opt tc=tc par=par;

proc iml;
use tetra;
read all into tetra;
close tetra;

start f(th) global(theta,tetra);
   theta = th;
   vv    = v(tetra[,1]);
   error = ssq(vv-tetra[,2]);
   return(error);
finish;

start v(t) global(theta);
   v = theta[3]*theta[1]/(theta[2]-theta[1])*
           (exp(-theta[1]*t)-exp(-theta[2]*t));
   return(v);
finish;

h      = {1.e-14 1. 1.e-5};
opt    = {0 2 0 1 };
tc     = repeat(.,1,12);
tc[1]  = 100;
tc[7]  = 1.e-8;
par    = { 1.e-13 . 1.e-10 . . . . };
con    = j(1,3,0.);
itheta = { .1  .3  10};
eps    = 1.e-11;

ods select ParameterEstimates#2;
call nlpdd(rc,rx,"f",itheta) blc=con opt=opt tc=tc par=par;

data oil(drop=temp);
   input temp time bitumen oil;
   cards;
673     5      0.      0.
673     7      2.2     0.
673    10     11.5      .7
673    15     13.7     7.2
673    20     15.1    11.5
673    25     17.3    15.8
673    30     17.3    20.9
673    40     20.1    26.6
673    50     20.1    32.4
673    60     22.3    38.1
673    80     20.9    43.2
673   100     11.5    49.6
673   120      6.5    51.8
673   150      3.6    54.7
;

proc iml;
use oil;
read all into a;
close oil;


/****************************************************************/
/* The INS function inserts a value given by FROM into a vector */
/* given by INTO, sorts the result, and posts the global matrix */
/* that can be used to delete the effects of the point FROM.    */
/****************************************************************/
start ins(from,into) global(permm);
   in    =  into // from;
   x     =  i(nrow(in));
   permm = inv(x[rank(in),]);
   return(permm*in);
finish;

start der(t,x) global(thmtrx,thet);
   y     = thmtrx*x;
   if ( t <= thet[5] )  then y = 0*y;
   return(y);
finish;

start jac(t,x) global(thmtrx,thet);
   y     = thmtrx;
   if ( t <= thet[5] )  then y = 0*y;
   return(y);
finish;

start f(theta) global(thmtrx,thet,time,h,a,eps,permm);
   thet = theta;
   thmtrx = (-(theta[1]+theta[4]) ||         0            || 0 )//
            (theta[1]             || -(theta[2]+theta[3]) || 0 )//
            (theta[4]             || theta[2]             || 0 );
   t = ins( theta[5],time);
   c = { 100, 0, 0};
   call ode( r1, "der",c , t , h) j="jac" eps=eps;

   /* send the intermediate value to the last column */
   r = (c ||r1) * permm;
   m = r[2:3,(2:nrow(time))];
   mm = m`- a[,2:3];
   call qr(q,r,piv,lindep,mm);
   v = det(r);
   return(abs(v));
finish;

opt = {0 2 0 1 };
tc = repeat(.,1,12);
tc[1] = 100;
tc[7] = 1.e-7;
par = { 1.e-13 . 1.e-10 . . . .};
con = j(1,5,0.);
h = {1.e-14 1. 1.e-5};
time = (0 // a[,1]);
eps = 1.e-5;
itheta = { 1.e-3 1.e-3 1.e-3 1.e-3 1.};

ods select ParameterEstimates#2;
call nlpqn(rc,rx,"f",itheta)  blc=con opt=opt tc=tc par=par;

proc iml;
use oil;
read all into a;
close oil;

start ins(from,into) global(permm);
   in    =  into // from;
   x     =  i(nrow(in));
   permm = inv(x[rank(in),]);
   return(permm*in);
finish;

start f(thet) global(time,a);
   do i = 1 to nrow(time);
      t   = time[i];
      v1  = 100;
      if ( t >= thet[5] ) then
         v1 = 100*ev(t,thet[1],thet[4],thet[5]);
      v2 = 0;
      if ( t >= thet[5] ) then
         v2 = 100*thet[1]/(thet[2]+thet[3]-thet[1]-thet[4])*
               (ev(t,thet[1],thet[4],thet[5])-
                ev(t,thet[2],thet[3],thet[5]));
      v3 = 0;
      if ( t >= thet[5] ) then
         v3  = 100*thet[4]/(thet[1]+thet[4])*
           (1. - ev(t,thet[1],thet[4],thet[5])) +
           100*thet[1]*thet[2]/(thet[2]+thet[3]-thet[1]-thet[4])*(
           (1.-ev(t,thet[1],thet[4],thet[5]))/(thet[1]+thet[4]) -
           (1.-ev(t,thet[2],thet[3],thet[5]))/(thet[2]+thet[3])    );
      y = y // (v1 || v2 || v3);
   end;
   mm = y[,2:3]-a[,2:3];
   call qr(q,r,piv,lindep,mm);
   v = det(r);
   return(abs(v));
finish;

start ev(t,a,b,c);
   return(exp(-(a+b)*(t-c)));
finish;

opt = {0 2 0 1 };
tc = repeat(.,1,12);
tc[1] = 100;
tc[7] = 1.e-7;
con     = { 0.  0.  0.  0.   .   .  . ,
             .   .   .   .   .    . . ,
            -1   1   1  -1   .   1  1.e-7 };
time    =  a[,1];
par     = { 1.e-13 . 1.e-10 . . . .};
itheta  = { 1.e-3 1.e-3 1.e-2 1.e-3 1.};

call nlpqn(rc,rx,"f",itheta)  blc=con opt=opt tc=tc par=par;