MODULE PmetricsFunctionModel%s;

	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 = %s;

	VAR
		equations: Equations;
		solver: MathODE.Solver;
		
			
	PROCEDURE InterpMoreCov3*(IN parameters: ARRAY OF REAL; t: REAL; firstindex: INTEGER; OUT p: ARRAY OF REAL);
		VAR
		n_cov_cont,  i,  flag, index, acc, endm, max_m, m, j: INTEGER;
		gradient: REAL;
		BEGIN
    	
		n_cov_cont:=MathFunc.Round(parameters[0]);
		
		FOR i:=1 TO n_cov_cont DO
			max_m:=MathFunc.Round(parameters[i]);
       		m:= MathFunc.Round(parameters[n_cov_cont+i]); 
			acc:=firstindex;
		
			IF (i#1) THEN
				FOR j:=1 TO (i-1) DO 
					acc:=acc+MathFunc.Round(parameters[i-j])*2;
				END;
			END;
		
	   		index:=acc;	
   			flag:=0;
  	 		WHILE (flag=0) & (index<m+acc-1) DO
      			IF (t >= parameters[index]) & (t <= parameters[index+1])THEN
        			endm:=index+1;
        			flag:=1;
      			ELSE
        			index:=index+1;
     	 		END;
			END;
		
   
      		gradient:= (t- parameters[endm-1])/(parameters[endm]- parameters[endm-1]);
     		p[i-1]:= parameters[max_m+endm-1]+ gradient*(parameters[max_m+endm] - parameters[max_m+endm-1]);
		END;
	END InterpMoreCov3;
		
			PROCEDURE InterpMoreCov3_CAT*(IN parameters: ARRAY OF REAL; t: REAL; firstindex: INTEGER; firstindex_cov: INTEGER; OUT p: ARRAY OF REAL); 
	VAR 
		n_cov_cat, i, flag, index, acc, endm, max_m, m, j: INTEGER; 
		gradient: REAL; 
	BEGIN 

		n_cov_cat:=MathFunc.Round(parameters[firstindex]);
		FOR i:=1 TO n_cov_cat DO 
			max_m:=MathFunc.Round(parameters[firstindex+i]); 
			m:= MathFunc.Round(parameters[firstindex+n_cov_cat+i]); 
			
			acc:=firstindex_cov; 
			IF (i#1) THEN 
				FOR j:=1 TO (i-1) DO 
					acc:=acc+MathFunc.Round(parameters[firstindex+i-j])*2; 
				END; 
			END; 

			index:=acc; 
			flag:=0; 
			WHILE (flag=0) & (index<=m+acc-2) DO 
				IF (t >= parameters[index]) & (t < parameters[index+1])THEN 
					p[i-1]:=parameters[max_m+index]; 
					flag:=1; 
				ELSIF (index=m+acc-2) & (t>=parameters[index+1]) THEN;
					p[i-1]:=parameters[max_m+index+1];
					flag:=1;
				ELSE
					index:=index+1; 
				END; 
			END;
		END; 
	END InterpMoreCov3_CAT;
		
				
	PROCEDURE UserDerivatives(IN %s, %s: ARRAY OF REAL;
	numEq: INTEGER; t: REAL; OUT d%sdt: ARRAY OF REAL) ;
	VAR
		n_cov_cont, n_cov_cat, lastindex, firstindex_cov_cont, firstindex_cov_cat, firstindex_cat, i: INTEGER; 
		%s 
		p_cont: ARRAY %s OF REAL; 
		p_cat: ARRAY %s OF REAL;
	BEGIN 
		n_cov_cont:=MathFunc.Round(theta[0]); 
		lastindex:=1+n_cov_cont*2;
		 
		firstindex_cat:=lastindex;
		n_cov_cat:=MathFunc.Round(theta[firstindex_cat]); 
		lastindex:=lastindex+1+n_cov_cat*2; 
		firstindex_cov_cont:=lastindex;
		
		FOR i:=1 TO n_cov_cont DO 
			lastindex:=MathFunc.Round(lastindex+2*MathFunc.Round(theta[i])); 
		END; 
		firstindex_cov_cat:=lastindex;

		FOR i:=1 TO n_cov_cat DO 
			lastindex:=MathFunc.Round(lastindex+2*MathFunc.Round(theta[i+n_cov_cont*2+1])); 
		END; 

		InterpMoreCov3(theta,t,firstindex_cov_cont, p_cont);
		%s

		InterpMoreCov3_CAT(theta,t,firstindex_cat, firstindex_cov_cat,p_cat);
		%s

		%s
		%s 
	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 := %s;
		m.F1Index := %s;
		m.tlag1Index := %s;
		m.nCmt := %s;
		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 PmetricsFunctionModel%s.
