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;