MODULE PmetricsFunctionModel1;

	IMPORT
		PmetricsPKModels,
		Math,
		MathODE,
		MathFunc,
		MathRungeKutta45,
		Vec := LibVectors,
		Solve := LibSolve;

	TYPE
		Equations = POINTER TO RECORD(MathODE.Equations) END;
		ODEModel* = POINTER TO RECORD (PmetricsPKModels.PKModel) END;
		SSFun  =  POINTER  TO  RECORD (Solve.RouteFns)
			m: PmetricsPKModels.PKModel; 
			p: PmetricsPKModels.ModelParameters;
			amt, rate, tau: REAL;
			cmt: INTEGER
		END;

	CONST tol = 1.0E-7;
		numEq = 6;

	VAR
		equations: Equations;
		solver: MathODE.Solver;

	PROCEDURE UserDerivatives(IN theta, DER99WB_UNIPV: ARRAY OF REAL;
	numEq: INTEGER; t: REAL; OUT dDER99WB_UNIPVdt: ARRAY OF REAL) ;
	VAR
		K1_UNIPV: REAL;
		K2_UNIPV: REAL;
		V1_UNIPV: REAL;
		C_UNIPV: REAL;
		PSI_UNIPV: REAL;
		K21_UNIPV: REAL;
		K12_UNIPV: REAL;
		LAMBDA0_UNIPV: REAL;
		LAMBDA1_UNIPV: REAL;
		K10_UNIPV: REAL;
		WTOT_UNIPV: REAL;
		K21_POP: REAL;
		K10_POP: REAL;
		K12_POP: REAL;
		PSI_POP: REAL;
		logLAMBDA0_POP: REAL;
		logLAMBDA1_POP: REAL;
		logK2_POP: REAL;
		V1_POP: REAL;
		logK1_POP: REAL;
		LAMBDA0_POP_UNIPV: REAL;
		LAMBDA1_POP_UNIPV: REAL;
		K2_POP_UNIPV: REAL;
		K1_POP_UNIPV: REAL;
		VAR0, VAR1, VAR2, VAR3, VAR4, VAR5: REAL;
		IC0, IC1, IC2, IC3, IC4, IC5: REAL;
		
	BEGIN
		K21_POP := theta[0];
		K10_POP := theta[1];
		K12_POP := theta[2];
		PSI_POP := theta[3];
		logLAMBDA0_POP := theta[4];
		logLAMBDA1_POP := theta[5];
		logK2_POP := theta[6];
		V1_POP := theta[7];
		logK1_POP := theta[8];
		IC0 := theta[9];
		IC1 := theta[10];
		IC2 := theta[11];
		IC3 := theta[12];
		IC4 := theta[13];
		IC5 := theta[14];
		

		VAR0 := DER99WB_UNIPV[0]+IC0;
		VAR1 := DER99WB_UNIPV[1]+IC1;
		VAR2 := DER99WB_UNIPV[2]+IC2;
		VAR3 := DER99WB_UNIPV[3]+IC3;
		VAR4 := DER99WB_UNIPV[4]+IC4;
		VAR5 := DER99WB_UNIPV[5]+IC5;
		
		K2_POP_UNIPV := (Math.Exp(logK2_POP));
		V1_UNIPV := V1_POP;
		K1_POP_UNIPV := (Math.Exp(logK1_POP));
		LAMBDA0_POP_UNIPV := (Math.Exp(logLAMBDA0_POP));
		PSI_UNIPV := PSI_POP;
		LAMBDA1_POP_UNIPV := (Math.Exp(logLAMBDA1_POP));
		K21_UNIPV := K21_POP;
		K12_UNIPV := K12_POP;
		LAMBDA0_UNIPV := LAMBDA0_POP_UNIPV;
		LAMBDA1_UNIPV := LAMBDA1_POP_UNIPV;
		WTOT_UNIPV := (((VAR2 + VAR3) + VAR4) + VAR5);
		K10_UNIPV := K10_POP;
		K1_UNIPV := K1_POP_UNIPV;
		K2_UNIPV := K2_POP_UNIPV;
		C_UNIPV := (VAR0 / V1_UNIPV);
		
		dDER99WB_UNIPVdt[0] := ((K21_UNIPV * VAR1) - ((K10_UNIPV + K12_UNIPV) * VAR0));
		dDER99WB_UNIPVdt[1] := ((K12_UNIPV * VAR0) - (K21_UNIPV * VAR1));
		dDER99WB_UNIPVdt[2] := (((LAMBDA0_UNIPV * VAR2) / MathFunc.Power((1 + MathFunc.Power(((WTOT_UNIPV * LAMBDA0_UNIPV) / LAMBDA1_UNIPV),PSI_UNIPV)),(1 / PSI_UNIPV))) - ((K2_UNIPV * C_UNIPV) * VAR2));
		dDER99WB_UNIPVdt[3] := (((K2_UNIPV * C_UNIPV) * VAR2) - (K1_UNIPV * VAR3));
		dDER99WB_UNIPVdt[4] := ((K1_UNIPV * VAR3) - (K1_UNIPV * VAR4));
		dDER99WB_UNIPVdt[5] := ((K1_UNIPV * VAR4) - (K1_UNIPV * VAR5));
			
	END UserDerivatives;

	PROCEDURE (equations: Equations) Derivatives* (IN theta, x: ARRAY OF REAL;
	numEq: INTEGER; t: REAL; OUT dxdt: ARRAY OF REAL) ;
	VAR
		i, iRate: INTEGER;
	BEGIN
		UserDerivatives (theta, x, numEq, t, dxdt);
		iRate := LEN(theta) - numEq;
		FOR i := 0 TO numEq-1 DO
			dxdt[i] := dxdt[i] + theta[iRate+i]
		END;
	
	END Derivatives;

	PROCEDURE (equations: Equations) SecondDerivatives* (IN theta, x: ARRAY OF REAL;
		numEq: INTEGER; t: REAL; OUT d2xdt2: ARRAY OF REAL);
	BEGIN
	END SecondDerivatives;

	PROCEDURE (equations: Equations) Jacobian* (IN theta, x: ARRAY OF REAL;
 	numEq: INTEGER; t: REAL; OUT jacob: ARRAY OF ARRAY OF REAL);
	BEGIN
	END Jacobian;
	
	PROCEDURE (m: ODEModel) InitModel*;
	BEGIN
		m.nParameter := 27;
		m.F1Index := 15;
		m.tlag1Index := 21;
		m.nCmt := 6;
		MathRungeKutta45.Install;
		solver := MathRungeKutta45.fact.New();
		NEW(equations);
		solver.Init(equations, numEq);
	END InitModel;

	PROCEDURE (m: ODEModel) Pred1*(dt: REAL; IN parameter: PmetricsPKModels.ModelParameters;
		IN init, rate: ARRAY OF REAL; OUT pred: ARRAY OF REAL);
	(* General nonlinear compartment model *)
	(* Calculate amount in each compartment after dt time units where init = initial conditions and 
	rate = rates into each compartment *) 
	VAR
		theta: POINTER TO ARRAY OF REAL;
		i, nCmt, np: INTEGER;
	BEGIN
		np := m.nParameter;
		nCmt := m.nCmt;
		NEW(theta,np+nCmt);
		FOR i := 0 TO np-1 DO
			theta[i] := parameter.parameters[i]
		END;
		FOR i := 0 TO nCmt-1 DO
			theta[np+i] := rate[i]
		END;
		
		solver.AccurateStep (theta, init, numEq, 0.0, dt, tol, pred);
	
	END Pred1;

	PROCEDURE (fn : SSFun) Eval (x : Vec.Vector) : Vec.Vector;
	VAR
		i, np: INTEGER;
		theta: POINTER TO ARRAY OF REAL;
		rate, xout, init: ARRAY PmetricsPKModels.nCmtMax OF REAL;
		result: Vec.Vector;
	BEGIN
		FOR i := 0 TO fn.m.nCmt-1 DO
			rate[i] := 0;
			init[i] := Math.Exp(x[i])
		END;
		
		result := Vec.New(fn.m.nCmt);
		IF fn.rate = 0 THEN	(* bolus dose *)
			init[fn.cmt-1] := init[fn.cmt-1] + fn.amt;
			fn.m.Pred1(fn.tau,fn.p,init,rate,xout);
			FOR i := 0 TO fn.m.nCmt-1 DO
				result[i] := xout[i] - Math.Exp(x[i])
			END;
		ELSIF fn.tau > 0 THEN	(* multiple truncated infusions *)
			rate[fn.cmt-1] := fn.rate;
			fn.m.Pred1(fn.amt/fn.rate,fn.p,init,rate,xout);
			rate[fn.cmt-1] := 0;
			FOR i := 0 TO fn.m.nCmt-1 DO
				init[i] := xout[i]
			END;
			fn.m.Pred1(fn.tau-fn.amt/fn.rate,fn.p,init,rate,xout);
			FOR i := 0 TO fn.m.nCmt-1 DO
				result[i] := xout[i] - Math.Exp(x[i])
			END;
		ELSE	(* constant infusion *)
			rate[fn.cmt-1] := fn.rate;
			np := fn.p.nParameter;
			NEW(theta,np);
			FOR i := 0 TO np-1 DO
				theta[i] := fn.p.parameters[i]
			END;
			UserDerivatives(theta, init, fn.m.nCmt, 0, xout);
			FOR i := 0 TO fn.m.nCmt-1 DO
				result[i] := xout[i] + rate[i]
			END;
		END;

		RETURN result;

	END  Eval;

	PROCEDURE (m:ODEModel) PredSS*(IN parameter: PmetricsPKModels.ModelParameters;
		amt, rate, ii: REAL; cmt: INTEGER; OUT pred: ARRAY OF REAL);
	(* General nonlinear compartment model *)
	(* Calculate amount in each compartment at the end of a steady-state dosing interval or during 
	a steady-state constant input (if ii = 0) *)
	
	(* Need to implement numerical solution to boundary value problem *)
	(* Find x0 such that pred = x0 - amt *)
	(* See LibSolve for root finding functions *)
	VAR
    	r, xBeg  :  Vec.Vector;
    	ssFun   :  SSFun;
    	route    :  Solve.Route;
		i, j, nCmt, np: INTEGER;
		zrate, init, xout: ARRAY PmetricsPKModels.nCmtMax OF REAL;
	BEGIN
		nCmt := m.nCmt;
		r := Vec.New (nCmt);
		FOR i := 0 TO nCmt-1 DO
			r[i] := 0
		END;
		xBeg   :=  Vec.New (nCmt);
		FOR i := 0 TO m.nCmt-1 DO
			zrate[i] := 0;
			init[i] := 0
		END;

		NEW (ssFun);
		ssFun.p := parameter;
		ssFun.amt := amt;
		ssFun.tau := ii;
		ssFun.cmt := cmt;
		ssFun.m := m;

		IF rate = 0 THEN	(* bolus dose *)
			
			init[cmt-1] := init[cmt-1] + amt;
			m.Pred1(ii,parameter,init,zrate,xout);
			FOR i := 0 TO nCmt-1 DO
				xBeg[i] := Math.Ln(xout[i])
			END;
			ssFun.rate := 0;
			route  :=  Solve.NewRoute (ssFun, xBeg, r, 0., 0., 0., 0., NIL);
			REPEAT  UNTIL  route.Leg ();
			WHILE ~(route.valTiny OR route.delTiny) DO
			(* This is a brute force approach to improving the initial estimates if the root *)
			(* solver does not converge *)
				FOR i := 1 TO 5 DO
					FOR j := 0 TO nCmt-1 DO
						init[j] := xout[j]
					END;
					init[cmt-1] := init[cmt-1] + amt;
					m.Pred1(ii,parameter,init,zrate,xout);
				END;
				FOR i := 0 TO nCmt-1 DO
					xBeg[i] := Math.Ln(xout[i])
				END;
				ssFun.rate := 0;
				route  :=  Solve.NewRoute (ssFun, xBeg, r, 0., 0., 0., 0., NIL);
				REPEAT  UNTIL  route.Leg ();
			END;
			ASSERT(route.valTiny OR route.delTiny,20);
			FOR i := 0 TO nCmt-1 DO
				pred[i] := Math.Exp(route.x[i])
			END;			

		ELSIF ii > 0 THEN	(* multiple truncated infusions *)

			zrate[cmt-1] := zrate[cmt-1] + rate;
			m.Pred1(amt/rate,parameter,init,zrate,xout);
			FOR i := 0 TO nCmt-1 DO
				init[i] := xout[i]
			END;
			zrate[cmt-1] := 0;
			m.Pred1(ii-amt/rate,parameter,init,zrate,xout);
			FOR i := 0 TO nCmt-1 DO
				xBeg[i] := Math.Ln(xout[i])
			END;
			ssFun.rate := rate;
			route  :=  Solve.NewRoute (ssFun, xBeg, r, 0., 0., 0., 0., NIL);
			REPEAT  UNTIL  route.Leg ();
			WHILE ~(route.valTiny OR route.delTiny) DO
			(* This is a brute force approach to improving the initial estimates if the root *)
			(* solver does not converge *)
				FOR i := 1 TO 5 DO
					FOR j := 0 TO nCmt-1 DO
						init[j] := xout[j]
					END;
					zrate[cmt-1] := zrate[cmt-1] + rate;
					m.Pred1(amt/rate,parameter,init,zrate,xout);
					zrate[cmt-1] := 0;
					m.Pred1(ii-amt/rate,parameter,init,zrate,xout);
				END;
				FOR i := 0 TO nCmt-1 DO
					xBeg[i] := Math.Ln(xout[i])
				END;
				ssFun.rate := rate;
				route  :=  Solve.NewRoute (ssFun, xBeg, r, 0., 0., 0., 0., NIL);
				REPEAT  UNTIL  route.Leg ();
			END;
			ASSERT(route.valTiny OR route.delTiny,20);
			FOR i := 0 TO nCmt-1 DO
				pred[i] := Math.Exp(route.x[i])
			END;
			
		ELSE	(* constant infusion *)

			zrate[cmt-1] := zrate[cmt-1] + rate;
			m.Pred1(100.0,parameter,init,zrate,xout);
			FOR i := 0 TO nCmt-1 DO
				xBeg[i] := Math.Ln(xout[i])
			END;
			ssFun.rate := rate;
			route  :=  Solve.NewRoute (ssFun, xBeg, r, 0., 0., 0., 0., NIL);
			REPEAT  UNTIL  route.Leg ();
			ASSERT(route.valTiny OR route.delTiny,20);
			FOR i := 0 TO nCmt-1 DO
				pred[i] := Math.Exp(route.x[i])
			END;
			
		END;
	
	END PredSS;

END PmetricsFunctionModel1.


