MODULE PmetricsPKModels;

	IMPORT Math, Strings, StdLog;

	CONST
		missing = -INF; nEventMax = %s; nParameterMax* = %s; nCmtMax* = %s;

	TYPE
		PKModel* = POINTER TO ABSTRACT RECORD
			nParameter*, F1Index*, tlag1Index*, nCmt*: INTEGER
		END;
		TwoCptModel* = POINTER TO RECORD (PKModel) END;
		OneCptModel* = POINTER TO RECORD (PKModel) END;
		
		String10* = ARRAY 10 OF CHAR;
		RealParameters* = ARRAY nParameterMax OF REAL;
		RealRates* = ARRAY nCmtMax OF REAL;
		EventHistory* = RECORD
			time, amt, rate, ii: ARRAY nEventMax OF REAL;
			evid, cmt, addl, ss: ARRAY nEventMax OF INTEGER;
			keep, new: ARRAY nEventMax OF BOOLEAN;
			nEvent*: INTEGER
		END;
		Event* = RECORD
			time*, amt*, rate*, ii*: REAL;
			evid*, cmt*, addl*, ss*: INTEGER;
			keep*, new*: BOOLEAN
		END;
		ModelParameterHistory* = RECORD
			time: ARRAY nEventMax OF REAL;
			parameters: ARRAY nEventMax OF RealParameters;
			nParameter: INTEGER;
			nEvent: INTEGER
		END;
		ModelParameters* = RECORD
			time*: REAL;
			parameters*: RealParameters;
			nParameter*: INTEGER
		END;
		RateHistory* = RECORD
			time: ARRAY nEventMax OF REAL;
			rate: ARRAY nEventMax OF RealRates;
			nEvent: INTEGER
		END;
		Rate* = RECORD
			time*: REAL;
			rate*: RealRates
		END;

	PROCEDURE  ClearEventHistory* (VAR events: EventHistory);
	BEGIN
		events.nEvent := 0;
	END ClearEventHistory;

	PROCEDURE  ClearRateHistory* (VAR rates: RateHistory);
	BEGIN
		rates.nEvent := 0;
	END ClearRateHistory;

	PROCEDURE  CheckEventHistory* (IN events: EventHistory): BOOLEAN;
	(* Checks whether EventHistory is in chronological order *)
	VAR
		i: INTEGER;
		ordered: BOOLEAN;
	BEGIN
		i := 1;
		ordered := TRUE;
		WHILE (i < events.nEvent-1) & (ordered) DO
			ordered := (events.time[i] >= events.time[i-1]) OR
				(events.evid[i] = 3) OR (events.evid[i] = 4);
			INC(i)
		END;
		RETURN ordered			
	END CheckEventHistory;

	PROCEDURE  CheckRateHistory* (IN rates: EventHistory): BOOLEAN;
	(* Checks whether EventHistory is in chronological order *)
	VAR
		i: INTEGER;
		ordered: BOOLEAN;
	BEGIN
		i := 0;
		ordered := TRUE;
		WHILE (i < rates.nEvent-1) & (ordered) DO
			ordered := rates.time[i] >= rates.time[i-1];
			INC(i)
		END;
		RETURN ordered			
	END CheckRateHistory;

	PROCEDURE GetEventI*(IN events: EventHistory; i: INTEGER; OUT event: Event);
	BEGIN
		event.time := events.time[i];
		event.evid := events.evid[i];
		event.cmt := events.cmt[i];
		event.amt := events.amt[i];
		event.rate := events.rate[i];
		event.addl := events.addl[i];
		event.ii := events.ii[i];
		event.ss := events.ss[i];
		event.keep := events.keep[i];
		event.new := events.new[i];
	END GetEventI;

	PROCEDURE GetEventValue*(IN events: EventHistory;
		i: INTEGER; name: ARRAY OF CHAR): REAL;
	VAR
		result: REAL;
	BEGIN
		IF name$ = "time" THEN
			result := events.time[i]
		ELSIF name$ = "evid" THEN
			result := events.evid[i]
		ELSIF name$ = "cmt" THEN
			result := events.cmt[i]
		ELSIF name$ = "amt" THEN
			result := events.amt[i]
		ELSIF name$ = "rate" THEN
			result := events.rate[i]
		ELSIF name$ = "addl" THEN
			result := events.addl[i]
		ELSIF name$ = "ii" THEN
			result := events.ii[i]
		ELSIF name$ = "ss" THEN
			result := events.ss[i]
		ELSIF name$ = "keep" THEN
			IF events.keep [i] THEN
				result := 1.0
			ELSE
				result := 0.0
			END
		ELSIF name$ = "new" THEN
			IF events.new [i] THEN
				result := 1.0
			ELSE
				result := 0.0
			END
		ELSE
			result := -INF
		END;
		RETURN result;
	END GetEventValue;

	PROCEDURE GetRateI*(IN rates: RateHistory; i: INTEGER; OUT rate: Rate);
	BEGIN
		rate.time := rates.time[i];
		rate.rate := rates.rate[i]
	END GetRateI;

	PROCEDURE GetModelParametersI*(IN parameters: ModelParameterHistory; i: INTEGER;
		OUT param: ModelParameters);
	BEGIN
		param.time := parameters.time[i];
		param.parameters := parameters.parameters[i];
		param.nParameter := parameters.nParameter
	END GetModelParametersI;
	
	PROCEDURE MakeEventHistory*(IN time, amt, rate, ii: ARRAY OF REAL;
		IN evid, cmt, addl, ss: ARRAY OF INTEGER; OUT events: EventHistory);
	VAR
		i: INTEGER;
	BEGIN
		(* length of time, amt, rate, evid and cmt must have same length *)
		(* ii, addl and ss must be length = 1 or the same length as time *)
		(* If ii, addl and ss have length = 1 then they are considered missing *) 
		events.nEvent := LEN(time);
		ASSERT((LEN(amt) = events.nEvent) & (LEN(rate) = events.nEvent) & (LEN(evid) = events.nEvent) & 
			(LEN(cmt) = events.nEvent), 20);
		ASSERT(((LEN(ii) = events.nEvent) OR (LEN(ii) = 1)) &
			((LEN(addl) = events.nEvent) OR (LEN(addl) = 1)) & 
			((LEN(ss) = events.nEvent) OR (LEN(ss) = 1)), 21);
		ASSERT(LEN(ii) = LEN(addl), 22);
		IF LEN(ii) = 1 THEN
			FOR i := 0 TO events.nEvent-1 DO
				events.ii[i] := 0;
				events.addl[i] := 0
			END;
		ELSE
			FOR i := 0 TO events.nEvent-1 DO
				events.ii[i] := ii[i];
				events.addl[i] := addl[i]
			END;
		END;
		IF LEN(ss) = 1 THEN
			FOR i := 0 TO events.nEvent-1 DO
				events.ss[i] := 0
			END;
		ELSE
			FOR i := 0 TO events.nEvent-1 DO
				events.ss[i] := ss[i]
			END;
		END;
		FOR i := 0 TO events.nEvent-1 DO
			events.time[i] := time[i];
			events.amt[i] := amt[i];
			events.rate[i] := rate[i];
			events.evid[i] := evid[i];
			events.cmt[i] := cmt[i];
			events.keep[i] := TRUE;
			events.new[i] := FALSE
		END;
	END MakeEventHistory;

	PROCEDURE MakeModelParameterHistory*(IN time, parameters: ARRAY OF REAL;
		nParameter: INTEGER; OUT parameterHistory: ModelParameterHistory);
	VAR
		i, j, np: INTEGER;
	BEGIN
		(* LEN(time) and LEN(parameters]/nParameter must be the same OR LEN(parameters]/nParameter = 1 *)
		np := LEN(parameters);
		ASSERT((np > 0) & (np MOD nParameter = 0), 21);
		parameterHistory.nEvent := np DIV nParameter;
		parameterHistory.nParameter := nParameter;
		ASSERT((LEN(time) = parameterHistory.nEvent) OR (parameterHistory.nEvent = 1), 21);
		FOR i := 0 TO parameterHistory.nEvent - 1 DO
			parameterHistory.time[i] := time[i];
			FOR j := 0 TO parameterHistory.nParameter - 1 DO
				parameterHistory.parameters[i,j] := parameters[nParameter*i + j];
			END;
		END;
	END MakeModelParameterHistory;

	PROCEDURE  InsertEventAtI* (VAR events: EventHistory; i: INTEGER; IN event: Event);
	VAR
		j: INTEGER;
	BEGIN
		ASSERT(0 <= i, 20);
		ASSERT(events.nEvent < nEventMax, 21);
		IF i > events.nEvent THEN
			i := events.nEvent
		END;
		FOR j := events.nEvent-1 TO i BY -1 DO
			events.time[j+1] := events.time[j];
			events.evid[j+1] := events.evid[j];
			events.cmt[j+1] := events.cmt[j];
			events.amt[j+1] := events.amt[j];
			events.rate[j+1] := events.rate[j];
			events.addl[j+1] := events.addl[j];
			events.ii[j+1] := events.ii[j];
			events.ss[j+1] := events.ss[j];
			events.keep[j+1] := events.keep[j];
			events.new[j+1] := events.new[j]
		END;
		events.time[i] := event.time;
		events.evid[i] := event.evid;
		events.cmt[i] := event.cmt;
		events.amt[i] := event.amt;
		events.rate[i] := event.rate;
		events.addl[i] := event.addl;
		events.ii[i] := event.ii;
		events.ss[i] := event.ss;
		events.keep[i] := event.keep;
		events.new[i] := event.new;
		INC(events.nEvent);
	END InsertEventAtI;

	PROCEDURE  InsertRateAtI* (VAR rates: RateHistory; i: INTEGER; IN rate: Rate);
	VAR
		j: INTEGER;
	BEGIN
		ASSERT(0 <= i, 20);
		ASSERT(rates.nEvent < nEventMax, 21);
		IF i > rates.nEvent THEN
			i := rates.nEvent
		END;
		FOR j := rates.nEvent-1 TO i BY -1 DO
			rates.time[j+1] := rates.time[j];
			rates.rate[j+1] := rates.rate[j];
		END;
		rates.time[i] := rate.time;
		rates.rate[i] := rate.rate;
		INC(rates.nEvent);
	END InsertRateAtI;

	PROCEDURE  InsertModelParametersAtI* (VAR parameters: ModelParameterHistory; i: INTEGER; 
		IN param: ModelParameters);
	VAR
		j: INTEGER;
	BEGIN
		ASSERT(0 <= i, 20);
		ASSERT(parameters.nEvent < nEventMax, 21);
		(* ASSERT(parameters.names = param.names, 22);  Can't do that *)
		IF i > parameters.nEvent THEN
			i := parameters.nEvent
		END;
		FOR j := parameters.nEvent-1 TO i BY -1 DO
			parameters.time[j+1] := parameters.time[j];
			parameters.parameters[j+1] := parameters.parameters[j];
		END;
		parameters.time[i] := param.time;
		parameters.parameters[i] := param.parameters;
		INC(parameters.nEvent);
	END InsertModelParametersAtI;

	PROCEDURE  RemoveEventI* (VAR events: EventHistory; i: INTEGER);
	VAR
		j: INTEGER;
	BEGIN
		ASSERT(0 <= i, 20);
		IF i <= events.nEvent-1 THEN
			FOR j := i TO events.nEvent-2 DO
				events.time[j] := events.time[j+1];
				events.evid[j] := events.evid[j+1];
				events.cmt[j] := events.cmt[j+1];	
				events.amt[j] := events.amt[j+1];
				events.rate[j] := events.rate[j+1];
				events.addl[j] := events.addl[j+1];
				events.ii[j] := events.ii[j+1];
				events.ss[j] := events.ss[j+1];
				events.keep[j] := events.keep[j+1];
				events.new[j] := events.new[j+1]
			END;
		END;
		DEC(events.nEvent);
	END RemoveEventI;
	
	PROCEDURE  RemoveRateI* (VAR rates: EventHistory; i: INTEGER);
	VAR
		j: INTEGER;
	BEGIN
		ASSERT(0 <= i, 20);
		IF i <= rates.nEvent-1 THEN
			FOR j := i TO rates.nEvent-2 DO
				rates.time[j] := rates.time[j+1];
				rates.rate[j] := rates.rate[j+1];
			END;
		END;
		DEC(rates.nEvent);
	END RemoveRateI;
	
	PROCEDURE CleanEvents*(VAR events: EventHistory);
	VAR
		i: INTEGER;
	BEGIN
		ASSERT(0 <= i, 20);
		FOR i := events.nEvent-1 TO 0 DO
			IF ~events.keep[i] THEN
				RemoveEventI(events, i)
			END;
		END;
	END CleanEvents;

	PROCEDURE SearchReal* (IN v: ARRAY OF REAL; numItm: INTEGER; srchNum: REAL;
			OUT i: INTEGER);
	(* returns the position of the largest value <= srchNum *)
		VAR
			first, mid, last: INTEGER;
	BEGIN
		ASSERT((0 <= numItm) & (numItm <= LEN(v)), 20);
		first := 0;
		last := numItm - 1;
		IF srchNum >= v[last] THEN
			i := last;
		ELSE
			WHILE first <= last DO
				mid := (first + last) DIV 2;
				IF srchNum < v[mid] THEN
					last := mid - 1
				ELSIF srchNum > v[mid+1] THEN
					first := mid + 1
				ELSE
					WHILE (srchNum = v[mid+1]) DO
						INC(mid)
					END;
					i := mid;
					RETURN
				END
			END;
		END;
	END SearchReal;
	
	PROCEDURE AddLagTimes*(VAR events: EventHistory; VAR parameters: ModelParameterHistory;
		IN tlagIndexes, tlagCmts: ARRAY OF INTEGER);
	(* Implements absorption lag times by modifying the times of the dosing events *)
	VAR
		nlag, pos, i, j, k, error, ipar: INTEGER;
		nonzero: ARRAY nCmtMax OF BOOLEAN;
		newEvent: Event;
	BEGIN
		nlag := LEN(tlagIndexes);
		ASSERT(LEN(tlagCmts) = nlag, 20);
		IF nlag > 0 THEN
			ASSERT((parameters.nEvent = events.nEvent) OR (parameters.nEvent = 1), 21);
			FOR i := events.nEvent-1 TO 0 BY -1 DO
				IF (events.evid[i] = 1) OR (events.evid[i] = 4) THEN
					j := 0;
					WHILE (events.cmt[i] # tlagCmts[j]) & (j < nlag) DO
						INC(j)
					END;
					ipar := MIN(i,parameters.nEvent-1);
					IF (events.cmt[i] = tlagCmts[j]) & (parameters.parameters[ipar,tlagIndexes[j]] # 0.0) THEN
						(* Assumes non-negative lag times *)
						ASSERT(parameters.parameters[ipar,tlagIndexes[j]] > 0.0, 22);
						GetEventI(events,i,newEvent);
						newEvent.time := newEvent.time + parameters.parameters[ipar,tlagIndexes[j]];
						newEvent.keep := FALSE;
						newEvent.new := TRUE;
						SearchReal(events.time, events.nEvent, newEvent.time,k);
						InsertEventAtI (events, k+1, newEvent);
						(* The following statement changes events so that CleanEvents does not return an *)
						(* object identical to the original *)
						events.evid[i] := 2
					END;
				END;
			END;
		END;
	END AddLagTimes;

	PROCEDURE AddlDoseEventsOld*(VAR events: EventHistory);
	VAR
		addlEvent, newEvent: Event;
		i, j, k: INTEGER;
	BEGIN
		i := 0;
		WHILE i < events.nEvent DO
			IF ((events.evid[i] = 1) OR (events.evid[i] = 4)) & (events.addl[i] > 0) & (events.ii[i] > 0) THEN
				GetEventI(events,i,addlEvent);
				newEvent := addlEvent;
				newEvent.addl := 0;
				newEvent.ii := 0;
				newEvent.ss := 0;
				newEvent.keep := FALSE;
				newEvent.new := TRUE;
				FOR j := 1 TO addlEvent.addl DO
					newEvent.time := addlEvent.time + j*addlEvent.ii;
					SearchReal(events.time, events.nEvent, newEvent.time,k);
					InsertEventAtI (events, k+1, newEvent)
				END;
			END;
			INC(i)
		END;
	END AddlDoseEventsOld;

	PROCEDURE AddlDoseEvents*(VAR events: EventHistory);
	VAR
		addlEvent, newEvent: Event;
		i, j, k: INTEGER;
	BEGIN
		i := 0;
		WHILE i < events.nEvent DO
			IF ((events.evid[i] = 1) OR (events.evid[i] = 4)) & (events.addl[i] > 0) & (events.ii[i] > 0) THEN
				GetEventI(events,i,addlEvent);
				newEvent := addlEvent;
				newEvent.addl := 0;
				newEvent.ii := 0;
				newEvent.ss := 0;
				newEvent.keep := FALSE;
				newEvent.new := TRUE;
				FOR j := 1 TO addlEvent.addl DO
					newEvent.time := addlEvent.time + j*addlEvent.ii;
					SearchReal(events.time, events.nEvent, newEvent.time,k);
					InsertEventAtI (events, k+1, newEvent)
				END;
			END;
			INC(i)
		END;
	END AddlDoseEvents;

	PROCEDURE MakeRates*(VAR events: EventHistory; VAR rates: RateHistory);
	(* Make infusion rate history and augment event history with end-of-infusion times *)
	VAR
		newEvent: Event;
		newRate: Rate;
		i, k, iRate: INTEGER;
		endTime: REAL;
	BEGIN
		FOR i := 0 TO nCmtMax-1 DO
			newRate.rate[i] := 0
		END;
		newRate.time := 0;
		ClearRateHistory(rates);
		
		FOR i := 0 TO events.nEvent-1 DO;
			newRate.time := events.time[i];
			InsertRateAtI (rates, i, newRate);
		END;
		ASSERT(rates.nEvent = events.nEvent,20);
		
		i := 0;
		WHILE i < events.nEvent DO
			IF ((events.evid[i] = 1) OR (events.evid[i] = 4)) & (events.rate[i] > 0) & (events.amt[i] > 0) THEN
				endTime := events.time[i] + events.amt[i]/events.rate[i];
				SearchReal(events.time, events.nEvent, endTime,k);
(*
				StdLog.String("k = "); StdLog.Int(k); StdLog.Ln;
*)
				GetEventI(events,k,newEvent);
				newEvent.time := endTime;
				newEvent.rate := 0;
				newEvent.amt := 0;
				newEvent.cmt := events.cmt[i];
				newEvent.evid := 2;
				newEvent.addl := 0;
				newEvent.ii := 0;
				newEvent.ss := 0;
				newEvent.keep := FALSE;
				newEvent.new := TRUE;
				InsertEventAtI (events, k+1, newEvent);
				newRate.time := endTime;
				InsertRateAtI (rates, k+1, newRate);
				FOR iRate := i+1 TO k+1 DO
					rates.rate[iRate,events.cmt[i]-1] := rates.rate[iRate,events.cmt[i]-1] + events.rate[i];
(*
					StdLog.String("time = "); StdLog.Real(rates.time[iRate]); StdLog.String("    ");
					StdLog.String("rate = "); StdLog.Real(rates.rate[iRate,events.cmt[i]-1]); StdLog.Ln;
*)
				END;
			END;
			INC(i)
		END;
	END MakeRates;
		
	PROCEDURE CompleteModelParameterHistory*(VAR parameters: ModelParameterHistory; VAR events: EventHistory);
	(* If parameters contains only a single set then it is replicated by the number of events in events *)
	(* Otherwise new parameter vectors are added corresponding to those in events with new = TRUE *)
	VAR
	 	newParameters: ModelParameters;
		i, iEvent: INTEGER;
	BEGIN
		ASSERT(parameters.nEvent > 0, 20);
		ASSERT(parameters.nParameter > 0, 20);
		ASSERT(parameters.nEvent <= events.nEvent, 21);
		iEvent := 0;
		FOR i := 0 TO parameters.nEvent-1 DO
			WHILE events.new[iEvent] DO
				INC(iEvent)
			END;
			ASSERT(parameters.time[i] = events.time[iEvent] ,22)   (* Equating REAL's may be too risky? *);
			INC(iEvent)
		END; 
		IF parameters.nEvent = 1 THEN
			FOR i := 1 TO events.nEvent-1 DO
				parameters.parameters[i] := parameters.parameters[0];
				parameters.time[i] := events.time[i];
				events.new[i] := FALSE;
			END;
			parameters.nEvent := events.nEvent
		ELSE
			iEvent := events.nEvent - 1;
			FOR i := parameters.nEvent-1 TO 0 BY -1 DO
				WHILE events.new[iEvent] DO
					GetModelParametersI(parameters, i, newParameters);
					newParameters.time := events.time[iEvent];
					InsertModelParametersAtI(parameters, i+1, newParameters);
					events.new[iEvent] := FALSE;
					DEC(iEvent)
				END;
				DEC(iEvent)
			END;
		END;
	END CompleteModelParameterHistory;

	PROCEDURE SearchString*(IN stringList: ARRAY OF ARRAY OF CHAR; numItm: INTEGER;
		srchString: ARRAY OF CHAR): INTEGER;
	VAR
		i: INTEGER;
	BEGIN
		ASSERT((0 <= numItm) & (numItm <= LEN(stringList)), 20);
		i := 0;
		WHILE (i <  numItm-1) & (stringList[i]$ # srchString$) DO
			INC(i)
		END;
		IF stringList[i]$ # srchString$ THEN
			i := -1
		END;
		RETURN i;
	END SearchString;

	PROCEDURE GetModelParameterValue*(IN parameters: ModelParameterHistory;
		iEvent, iParameter: INTEGER): REAL;
		BEGIN
			ASSERT((iEvent >= 0) & (iEvent < parameters.nEvent), 20);
			ASSERT((iParameter >= 0) & (iParameter < parameters.nParameter) ,21);
			RETURN parameters.parameters[iEvent,iParameter]
		END GetModelParameterValue;

	PROCEDURE PolyExp*(x,dose,rate,xinf,tau: REAL; ss: BOOLEAN; a,alpha: ARRAY OF REAL;
		n: INTEGER): REAL;
	VAR
		result,bolusResult,dx,nIntv: REAL;
		i: INTEGER;
	BEGIN
		ASSERT((LEN(alpha) >= n) & (LEN(a) >= n), 20);
		result := 0;
		IF dose > 0 THEN (* bolus dose *)
			IF tau <=0 THEN
				IF x >= 0 THEN
					FOR i := 0 TO n-1 DO;
						result := result + a[i] * Math.Exp(-alpha[i]*x)
					END;
				END;
			ELSIF ~ss THEN
				nIntv := ENTIER(x/tau) + 1;
				dx := x - ENTIER(x/tau)*tau;
				FOR i := 0 TO n-1 DO;
					result := result + a[i] * Math.Exp(-alpha[i]*x)*
						(1-Math.Exp(-nIntv*alpha[i]*tau))/(1-Math.Exp(-alpha[i]*tau))
				END;
			ELSE
				dx := x - ENTIER(x/tau)*tau;
				FOR i := 0 TO n-1 DO;
					result := result + a[i] * Math.Exp(-alpha[i]*x)/(1-Math.Exp(-alpha[i]*tau))
				END;
			END;
		END;
		bolusResult := dose * result;
		result := 0;
		IF (rate > 0) & (xinf < INF) THEN (* truncated infusion *)
			IF tau <=0 THEN
				IF x >= 0 THEN
					IF x <= xinf THEN
						FOR i := 0 TO n-1 DO;
							result := result + a[i] * (1 - Math.Exp(-alpha[i]*x)) / alpha[i]
						END;
					ELSE
						FOR i := 0 TO n-1 DO;
							result := result + a[i] * (1 - Math.Exp(-alpha[i]*xinf))*Math.Exp(-alpha[i]*(x-xinf)) / alpha[i]
						END;
					END;
				END;
			ELSIF ~ss THEN
				ASSERT(xinf <= tau, 20); (* add other case later *)
				dx := x - ENTIER(x/tau)*tau;
				nIntv := ENTIER(x/tau) + 1;
				IF dx <= xinf THEN
					FOR i := 0 TO n-1 DO;
						IF n > 1 THEN
							result := result + a[i] * (1 - Math.Exp(-alpha[i]*xinf))*Math.Exp(-alpha[i]*(dx-xinf+tau)) *
								(1-Math.Exp(-(nIntv-1)*alpha[i]*tau))/(1-Math.Exp(-alpha[i]*tau)) / alpha[i]
						END;
						result := result + a[i] * (1 - Math.Exp(-alpha[i]*dx)) / alpha[i]
					END;
				ELSE
					FOR i := 0 TO n-1 DO;
						result := result + a[i] * (1 - Math.Exp(-alpha[i]*xinf))*Math.Exp(-alpha[i]*(dx-xinf)) *
							(1-Math.Exp(-nIntv*alpha[i]*tau))/(1-Math.Exp(-alpha[i]*tau)) / alpha[i]
					END;
				END;
			ELSE
				ASSERT(xinf <= tau, 20); (* add other case later *)
				dx := x - ENTIER(x/tau)*tau;
				nIntv := ENTIER(x/tau) + 1;
				IF dx <= xinf THEN
					FOR i := 0 TO n-1 DO;
						result := result + a[i] * (1 - Math.Exp(-alpha[i]*xinf))*Math.Exp(-alpha[i]*(dx-xinf+tau)) /
							(1-Math.Exp(-alpha[i]*tau)) / alpha[i] + a[i] * (1 - Math.Exp(-alpha[i]*dx)) / alpha[i]
					END;
				ELSE
					FOR i := 0 TO n-1 DO;
						result := result + a[i] * (1 - Math.Exp(-alpha[i]*xinf))*Math.Exp(-alpha[i]*(dx-xinf)) /
							(1-Math.Exp(-alpha[i]*tau)) / alpha[i]
					END;
				END;
			END;
		ELSE (* continuous infusion (xinf = INF). tau is ignored *)
			IF ~ss THEN
				IF x >= 0 THEN
					FOR i := 0 TO n-1 DO;
						result := result + a[i] * (1 - Math.Exp(-alpha[i]*x)) / alpha[i]
					END;
				END;
			ELSE
				FOR i := 0 TO n-1 DO;
					result := result + a[i]/alpha[i]
				END;
			END;
		END;
		RETURN bolusResult + rate * result;

	END PolyExp;

	PROCEDURE (m: PKModel) GetNParameter(): INTEGER, NEW;
	BEGIN
		RETURN m.nParameter;
	END GetNParameter;
	
	PROCEDURE (m: PKModel) GetF1Index(): INTEGER, NEW;
	BEGIN
		RETURN m.F1Index;
	END GetF1Index;
	
	PROCEDURE (m: PKModel) GetTlag1Index(): INTEGER, NEW;
	BEGIN
		RETURN m.tlag1Index;
	END GetTlag1Index;
	
	PROCEDURE (m: PKModel) GetNCmt(): INTEGER, NEW;
	BEGIN
		RETURN m.nCmt;
	END GetNCmt;
	
	PROCEDURE (m: PKModel) InitModel*, NEW, ABSTRACT;

	PROCEDURE (m: PKModel) Pred1*(dt: REAL; IN parameter: ModelParameters;
		IN init, rate: ARRAY OF REAL; OUT pred: ARRAY OF REAL), NEW, ABSTRACT;

	PROCEDURE (m: PKModel) PredSS*(IN parameter: ModelParameters; amt, rate, ii: REAL;
		cmt: INTEGER; OUT pred: ARRAY OF REAL), NEW, ABSTRACT;

	PROCEDURE (m: PKModel) Pred*(IN pArray,time,amt,rate,ii: ARRAY OF REAL;
		IN evid,cmt,addl,ss: ARRAY OF INTEGER; OUT pred: ARRAY OF REAL), NEW;
	(* Two compartment model with 1st order absorption *)
	(* Calculate amount in each compartment*) 

		VAR
			i, j, np, ikeep, nParameter, nCmt, F1Index, tlag1Index: INTEGER;
			dt, tprev: REAL;
			init: ARRAY nCmtMax OF REAL;
			tlagIndexes, tlagCmts: ARRAY nCmtMax OF INTEGER;
			pred1: ARRAY nCmtMax OF REAL;
			event: Event;
			rate2: Rate;
			events: EventHistory;
			parameters: ModelParameterHistory;
			parameter: ModelParameters;
			rates: RateHistory;

		BEGIN
			m.InitModel;
			nParameter := m.GetNParameter();
			nCmt := m.GetNCmt();
			F1Index := m.GetF1Index();
			tlag1Index := m.GetTlag1Index();
			MakeEventHistory(time,amt,rate,ii,evid,cmt,addl,ss,events);
			np := LEN(pArray);
			ASSERT((np > 0) & (np MOD nParameter = 0), 20);
			MakeModelParameterHistory(time,pArray,nParameter,parameters);
			ASSERT(CheckEventHistory(events), 21);
			AddlDoseEvents(events);
			CompleteModelParameterHistory(parameters, events);
			FOR i := 0 TO nCmt-1 DO
				tlagIndexes[i] := tlag1Index + i;
				tlagCmts[i] := i + 1;
			END;
			AddLagTimes(events, parameters, tlagIndexes, tlagCmts);
			MakeRates(events, rates);
(*			
			FOR i := 0 TO rates.nEvent-1 DO;
				StdLog.Real(rates.time[i]); StdLog.String("    "); StdLog.Real(rates.rate[i,0]); StdLog.Ln;
			END;
*)			
			CompleteModelParameterHistory(parameters, events);
			FOR i := 0 TO nCmt-1 DO
				init[i] := 0.0;
			END;
			GetEventI(events,0,event);
			tprev := event.time;
			ikeep := 0;
			FOR i := 0 TO events.nEvent-1 DO
				GetEventI(events,i,event);
				GetRateI(rates,i,rate2);
				FOR j := 0 TO nCmt-1 DO
(*
					StdLog.Int(j); StdLog.String("    rate2 ="); StdLog.Real(rate2.rate[j]); StdLog.Ln;
*)
					rate2.rate[j] := GetModelParameterValue(parameters, i, F1Index+j) * rate2.rate[j]
				END;
				GetModelParametersI(parameters, i, parameter);
				IF (event.evid = 3) OR (event.evid = 4) THEN	(* reset dose event *)
					dt := 0.0;
					FOR j := 0 TO nCmt-1 DO
						init[j] := 0.0;
					END;
				ELSE
					dt := event.time - tprev;
					m.Pred1(dt, parameter, init, rate2.rate, pred1);
					FOR j := 0 TO nCmt - 1 DO
						init[j] := pred1[j];
					END;
				END;
				IF ((event.evid = 1) OR (event.evid = 4)) &
					((event.ss = 1) OR (event.ss = 2) OR (event.ss = 3)) THEN	(* steady-state dose event *)
					m.PredSS(parameter, 
						GetModelParameterValue(parameters, i, F1Index+event.cmt-1) * event.amt,
						event.rate, event.ii, event.cmt, pred1); 
					IF event.ss = 2 THEN	(* steady-state without reset *)
						FOR j := 0 TO nCmt - 1 DO
							init[j] := init[j] + pred1[j];
						END;
					ELSE	(* steady-state with reset (ss = 1) *)
						FOR j := 0 TO nCmt - 1 DO
							init[j] := pred1[j];
						END;
					END;
				END;				
				IF ((event.evid = 1) OR (event.evid = 4)) & (event.rate = 0) THEN	(* bolus dose *)
					init[event.cmt-1] := init[event.cmt-1] + 
						GetModelParameterValue(parameters, i, F1Index+event.cmt-1) * event.amt
				END;
				IF event.keep THEN
					FOR j := 0 TO nCmt - 1 DO
						pred[nCmt*ikeep+j] := init[j];
(*						StdLog.Real(pred[nCmt*ikeep+j]); StdLog.Ln; *)
					END;
					INC(ikeep);
				END;
				tprev := event.time;
			END;

	END Pred;

	PROCEDURE (m: TwoCptModel) InitModel*;
	BEGIN
		m.nParameter := 11;
		m.F1Index := 5;
		m.tlag1Index := 8;
		m.nCmt := 3
	END InitModel;

	PROCEDURE (m: TwoCptModel) Pred1*(dt: REAL; IN parameter: ModelParameters;
		IN init, rate: ARRAY OF REAL; OUT pred: ARRAY OF REAL);
	(* Two compartment model with 1st order absorption *)
	(* Calculate amount in each compartment after dt time units where init = initial conditions and 
	rate = rates into each compartment *) 
	VAR
		a, alpha: ARRAY 3 OF REAL;
		CL, Q, V2, V3, dka, ka, k10, k12, k21, ksum: REAL;
		i: INTEGER;
	BEGIN
(*
		StdLog.Real(dt); StdLog.String("   "); StdLog.Real(init[0]); StdLog.String("   "); 
		StdLog.Real(rate[0]); .Ln;
*)
		CL := parameter.parameters[0];
		Q := parameter.parameters[1];
		V2 := parameter.parameters[2];
		V3 := parameter.parameters[3];
		dka := parameter.parameters[4];	
		ASSERT((CL > 0) & (Q > 0) & (V2 > 0) & (V3 > 0) & (dka > 0) , 20);
		k10 := CL/V2;
		k12 := Q/V2;
		k21 := Q/V3;
		ksum := k10 + k12 + k21;
		alpha[0] := (ksum + Math.Sqrt(ksum*ksum-4.0*k10*k21))/2.0;
		alpha[1] := (ksum - Math.Sqrt(ksum*ksum-4.0*k10*k21))/2.0;
		ka := dka + alpha[1];
		alpha[2] := ka;

		pred[0] := 0;
		pred[1] := 0;
		pred[2] := 0;

		IF (init[0] # 0.0) OR (rate[0] # 0.0) THEN
			pred[0] := init[0] * Math.Exp(-ka*dt);
			pred[0] := pred[0] + rate[0] * (1.0 - Math.Exp(-ka*dt)) / ka;
			a[0] := ka * (k21 - alpha[0])/((ka-alpha[0])*(alpha[1]-alpha[0])) ; 
			a[1] := ka * (k21 - alpha[1])/((ka-alpha[1])*(alpha[0]-alpha[1])) ; 
			a[2] := -(a[0] + a[1]);
			pred[1] := pred[1] +
				PolyExp(dt,init[0],0,0,0,FALSE,a,alpha,3) +
				PolyExp(dt,0,rate[0],dt,0,FALSE,a,alpha,3);
			a[0] := ka * k12/((ka-alpha[0])*(alpha[1]-alpha[0])) ; 
			a[1] := ka * k12/((ka-alpha[1])*(alpha[0]-alpha[1])) ; 
			a[2] := -(a[0] + a[1]);
			pred[2] := pred[2] +
				PolyExp(dt,init[0],0,0,0,FALSE,a,alpha,3) +
				PolyExp(dt,0,rate[0],dt,0,FALSE,a,alpha,3);
		END;
		
		IF (init[1] # 0.0) OR (rate[1] # 0.0) THEN
			a[0] := (k21 - alpha[0])/(alpha[1]-alpha[0]) ; 
			a[1] := (k21 - alpha[1])/(alpha[0]-alpha[1]) ; 
			pred[1] := pred[1] +
				PolyExp(dt,init[1],0,0,0,FALSE,a,alpha,2) +
				PolyExp(dt,0,rate[1],dt,0,FALSE,a,alpha,2);
			a[0] := k12/(alpha[1]-alpha[0]) ; 
			a[1] := -a[0]; 
			pred[2] := pred[2] +
				PolyExp(dt,init[1],0,0,0,FALSE,a,alpha,2) +
				PolyExp(dt,0,rate[1],dt,0,FALSE,a,alpha,2);
		END;

		IF (init[2] # 0.0) OR (rate[2] # 0.0) THEN
			a[0] := k21/(alpha[1]-alpha[0]) ; 
			a[1] := -a[0]; 
			pred[1] := pred[1] +
				PolyExp(dt,init[2],0,0,0,FALSE,a,alpha,2) +
				PolyExp(dt,0,rate[2],dt,0,FALSE,a,alpha,2);
			a[0] := (k10 + k12 - alpha[0])/(alpha[1]-alpha[0]) ; 
			a[1] := (k10 + k12 - alpha[1])/(alpha[0]-alpha[1]) ; 
			pred[2] := pred[2] +
				PolyExp(dt,init[2],0,0,0,FALSE,a,alpha,2) +
				PolyExp(dt,0,rate[2],dt,0,FALSE,a,alpha,2);
		END;
	
	END Pred1;

	PROCEDURE (m:TwoCptModel) PredSS*(IN parameter: ModelParameters; amt, rate, ii: REAL;
		cmt: INTEGER; OUT pred: ARRAY OF REAL);
	(* Two compartment model with 1st order absorption *)
	(* Calculate amount in each compartment at the end of a steady-state dosing interval or during 
	a steady-state constant input (if ii = 0) *)
	VAR
		a, alpha: ARRAY 3 OF REAL;
		CL, Q, V2, V3, dka, ka, k10, k12, k21, ksum: REAL;
		i: INTEGER;
	BEGIN
		CL := parameter.parameters[0];
		Q := parameter.parameters[1];
		V2 := parameter.parameters[2];
		V3 := parameter.parameters[3];
		dka := parameter.parameters[4];	
		ASSERT((CL > 0) & (Q > 0) & (V2 > 0) & (V3 > 0) & (dka > 0) , 20);
		k10 := CL/V2;
		k12 := Q/V2;
		k21 := Q/V3;
		ksum := k10 + k12 + k21;
		alpha[0] := (ksum + Math.Sqrt(ksum*ksum-4.0*k10*k21))/2.0;
		alpha[1] := (ksum - Math.Sqrt(ksum*ksum-4.0*k10*k21))/2.0;
		ka := dka + alpha[1];
		alpha[2] := ka;
		
		pred[0] := 0;
		pred[1] := 0;
		pred[2] := 0;

		IF rate = 0 THEN	(* bolus dose *)

			IF (cmt = 1) THEN
				a[0] := 0.0;
				a[1] := 0.0;
				a[2] := 1.0;
				pred[0] := PolyExp(ii,amt,0,0,ii,TRUE,a,alpha,3);
				a[0] := ka * (k21 - alpha[0])/((ka-alpha[0])*(alpha[1]-alpha[0])) ; 
				a[1] := ka * (k21 - alpha[1])/((ka-alpha[1])*(alpha[0]-alpha[1])) ; 
				a[2] := -(a[0] + a[1]);
				pred[1] := PolyExp(ii,amt,0,0,ii,TRUE,a,alpha,3);
				a[0] := ka * k12/((ka-alpha[0])*(alpha[1]-alpha[0])) ; 
				a[1] := ka * k12/((ka-alpha[1])*(alpha[0]-alpha[1])) ; 
				a[2] := -(a[0] + a[1]);
				pred[2] := PolyExp(ii,amt,0,0,ii,TRUE,a,alpha,3);
		
			ELSIF (cmt = 2) THEN
				a[0] := (k21 - alpha[0])/(alpha[1]-alpha[0]) ; 
				a[1] := (k21 - alpha[1])/(alpha[0]-alpha[1]) ; 
				pred[1] := PolyExp(ii,amt,0,0,ii,TRUE,a,alpha,2);
				a[0] := k12/(alpha[1]-alpha[0]) ; 
				a[1] := -a[0]; 
				pred[2] := PolyExp(ii,amt,0,0,ii,TRUE,a,alpha,2);

			ELSE	(* cmt = 3 *)
				a[0] := k21/(alpha[1]-alpha[0]) ; 
				a[1] := -a[0]; 
				pred[1] := PolyExp(ii,amt,0,0,ii,TRUE,a,alpha,2);
				a[0] := (k10 + k12 - alpha[0])/(alpha[1]-alpha[0]) ; 
				a[1] := (k10 + k12 - alpha[1])/(alpha[0]-alpha[1]) ; 
				pred[2] := PolyExp(ii,amt,0,0,ii,TRUE,a,alpha,2);
			END;
	
		ELSIF ii > 0 THEN	(* multiple truncated infusions *)

			IF (cmt = 1) THEN
				a[0] := 0.0;
				a[1] := 0.0;
				a[2] := 1.0;
				pred[0] := PolyExp(ii,0,rate,amt/rate,ii,TRUE,a,alpha,3);
				a[0] := ka * (k21 - alpha[0])/((ka-alpha[0])*(alpha[1]-alpha[0])) ; 
				a[1] := ka * (k21 - alpha[1])/((ka-alpha[1])*(alpha[0]-alpha[1])) ; 
				a[2] := -(a[0] + a[1]);
				pred[1] := PolyExp(ii,0,rate,amt/rate,ii,TRUE,a,alpha,3);
				a[0] := ka * k12/((ka-alpha[0])*(alpha[1]-alpha[0])) ; 
				a[1] := ka * k12/((ka-alpha[1])*(alpha[0]-alpha[1])) ; 
				a[2] := -(a[0] + a[1]);
				pred[2] := PolyExp(ii,0,rate,amt/rate,ii,TRUE,a,alpha,3);
		
			ELSIF (cmt = 2) THEN
				a[0] := (k21 - alpha[0])/(alpha[1]-alpha[0]) ; 
				a[1] := (k21 - alpha[1])/(alpha[0]-alpha[1]) ; 
				pred[1] := PolyExp(ii,0,rate,amt/rate,ii,TRUE,a,alpha,2);
				a[0] := k12/(alpha[1]-alpha[0]) ; 
				a[1] := -a[0]; 
				pred[2] := PolyExp(ii,0,rate,amt/rate,ii,TRUE,a,alpha,2);

			ELSE	(* cmt = 3 *)
				a[0] := k21/(alpha[1]-alpha[0]) ; 
				a[1] := -a[0]; 
				pred[1] := PolyExp(ii,0,rate,amt/rate,ii,TRUE,a,alpha,2);
				a[0] := (k10 + k12 - alpha[0])/(alpha[1]-alpha[0]) ; 
				a[1] := (k10 + k12 - alpha[1])/(alpha[0]-alpha[1]) ; 
				pred[2] := PolyExp(ii,0,rate,amt/rate,ii,TRUE,a,alpha,2);
			END;
			
		ELSE	(* constant infusion *)
		
			IF (cmt = 1) THEN
				a[0] := 0.0;
				a[1] := 0.0;
				a[2] := 1.0;
				pred[0] := PolyExp(0,0,rate,INF,0,TRUE,a,alpha,3);
				a[0] := ka * (k21 - alpha[0])/((ka-alpha[0])*(alpha[1]-alpha[0])) ; 
				a[1] := ka * (k21 - alpha[1])/((ka-alpha[1])*(alpha[0]-alpha[1])) ; 
				a[2] := -(a[0] + a[1]);
				pred[1] := PolyExp(0,0,rate,INF,0,TRUE,a,alpha,3);
				a[0] := ka * k12/((ka-alpha[0])*(alpha[1]-alpha[0])) ; 
				a[1] := ka * k12/((ka-alpha[1])*(alpha[0]-alpha[1])) ; 
				a[2] := -(a[0] + a[1]);
				pred[2] := PolyExp(0,0,rate,INF,0,TRUE,a,alpha,3);
		
			ELSIF (cmt = 2) THEN
				a[0] := (k21 - alpha[0])/(alpha[1]-alpha[0]) ; 
				a[1] := (k21 - alpha[1])/(alpha[0]-alpha[1]) ; 
				pred[1] := PolyExp(0,0,rate,INF,0,TRUE,a,alpha,2);
				a[0] := k12/(alpha[1]-alpha[0]) ; 
				a[1] := -a[0]; 
				pred[2] := PolyExp(0,0,rate,INF,0,TRUE,a,alpha,2);

			ELSE	(* cmt = 3 *)
				a[0] := k21/(alpha[1]-alpha[0]) ; 
				a[1] := -a[0]; 
				pred[1] := PolyExp(0,0,rate,INF,0,TRUE,a,alpha,2);
				a[0] := (k10 + k12 - alpha[0])/(alpha[1]-alpha[0]) ; 
				a[1] := (k10 + k12 - alpha[1])/(alpha[0]-alpha[1]) ; 
				pred[2] := PolyExp(0,0,rate,INF,0,TRUE,a,alpha,2);
			END;
		
		END;

	END PredSS;

(*****************************************************************************************************)
	
	PROCEDURE (m: OneCptModel) InitModel*;
	BEGIN
		m.nParameter := 7;
		m.F1Index := 3;
		m.tlag1Index := 5;
		m.nCmt := 2
	END InitModel;

	PROCEDURE (m: OneCptModel) Pred1*(dt: REAL; IN parameter: ModelParameters;
		IN init, rate: ARRAY OF REAL; OUT pred: ARRAY OF REAL);
	(* One compartment model with 1st order absorption *)
	(* Calculate amount in each compartment after dt time units where init = initial conditions and 
	rate = rates into each compartment *) 
	VAR
		a, alpha: ARRAY 2 OF REAL;
		CL, V2, dka, ka, k10: REAL;
		i: INTEGER;
	BEGIN
		CL := parameter.parameters[0];
		V2 := parameter.parameters[1];
		dka := parameter.parameters[2];	
		ASSERT((CL > 0) & (V2 > 0) & (dka > 0) , 20);
		k10 := CL/V2;
		alpha[0] := k10;
		ka := dka + alpha[0];
		alpha[1] := ka;

		pred[0] := 0;
		pred[1] := 0;

		IF (init[0] # 0.0) OR (rate[0] # 0.0) THEN
			pred[0] := init[0] * Math.Exp(-ka*dt);
			pred[0] := pred[0] + rate[0] * (1.0 - Math.Exp(-ka*dt)) / ka;
			a[0] := ka / (ka-alpha[0]) ; 
			a[1] := -a[0];
			pred[1] := pred[1] +
				PolyExp(dt,init[0],0,0,0,FALSE,a,alpha,2) +
				PolyExp(dt,0,rate[0],dt,0,FALSE,a,alpha,2);
		END;
		
		IF (init[1] # 0.0) OR (rate[1] # 0.0) THEN
			a[0] := 1.0 ; 
			pred[1] := pred[1] +
				PolyExp(dt,init[1],0,0,0,FALSE,a,alpha,1) +
				PolyExp(dt,0,rate[1],dt,0,FALSE,a,alpha,1);
		END;
	
	END Pred1;

	PROCEDURE (m:OneCptModel) PredSS*(IN parameter: ModelParameters; amt, rate, ii: REAL;
		cmt: INTEGER; OUT pred: ARRAY OF REAL);
	(* One compartment model with 1st order absorption *)
	(* Calculate amount in each compartment at the end of a steady-state dosing interval or during 
	a steady-state constant input (if ii = 0) *)
	VAR
		a, alpha: ARRAY 2 OF REAL;
		CL, V2, dka, ka, k10: REAL;
		i: INTEGER;
	BEGIN
		CL := parameter.parameters[0];
		V2 := parameter.parameters[2];
		dka := parameter.parameters[4];	
		ASSERT((CL > 0) & (V2 > 0) & (dka > 0) , 20);
		k10 := CL/V2;
		alpha[0] := k10;
		ka := dka + alpha[0];
		alpha[1] := ka;
		
		pred[0] := 0;
		pred[1] := 0;

		IF rate = 0 THEN	(* bolus dose *)

			IF (cmt = 1) THEN
				a[0] := 0.0;
				a[1] := 1.0;
				pred[0] := PolyExp(ii,amt,0,0,ii,TRUE,a,alpha,2);
				a[0] := ka / (ka-alpha[0]) ; 
				a[1] := -a[0];
				pred[1] := PolyExp(ii,amt,0,0,ii,TRUE,a,alpha,2);
		
			ELSE (* cmt = 2 *)
				a[0] := 1.0 ; 
				pred[1] := PolyExp(ii,amt,0,0,ii,TRUE,a,alpha,1);
			END;
	
		ELSIF ii > 0 THEN	(* multiple truncated infusions *)

			IF (cmt = 1) THEN
				a[0] := 0.0;
				a[1] := 1.0;
				pred[0] := PolyExp(ii,0,rate,amt/rate,ii,TRUE,a,alpha,2);
				a[0] := ka / (ka-alpha[0]) ; 
				a[1] := -a[0];
				pred[1] := PolyExp(ii,0,rate,amt/rate,ii,TRUE,a,alpha,2);
		
			ELSE (* cmt = 2 *)
				a[0] := 1.0 ; 
				pred[1] := PolyExp(ii,0,rate,amt/rate,ii,TRUE,a,alpha,1);
			END;
			
		ELSE	(* constant infusion *)
		
			IF (cmt = 1) THEN
				a[0] := 0.0;
				a[1] := 1.0;
				pred[0] := PolyExp(0,0,rate,0,0,TRUE,a,alpha,2);
				a[0] := ka / (ka-alpha[0]) ; 
				a[1] := -a[0];
				pred[1] := PolyExp(0,0,rate,0,0,TRUE,a,alpha,2);
		
			ELSE (* cmt = 2 *)
				a[0] := 1.0 ; 
				pred[1] := PolyExp(0,0,rate,0,0,TRUE,a,alpha,1);
			END;
		
		END;

	END PredSS;
		
	END PmetricsPKModels.
	
