.packageName <- "PK"


AUC <- function(time, conc, exact=NA, numintp=2, numtail=3, prev=0) {		     

	# function for linear interpolation/extrapolation
	linpol <- function(conc, time, exact){
		parms <- lm(conc~time)
		return(parms$coef[2] * exact + parms$coef[1])
	}
         
	# function to add parts of auc and aumc
	add <- function(time, conc) {
		auc <- 0; aumc <- 0
		for (i in 2:length(conc)) {
			auc <- auc  + 1/2 * (time[i]-time[i-1]) * (conc[i]+conc[i-1])
			aumc <- aumc + 1/2 * (time[i]-time[i-1]) * (conc[i]*time[i] + conc[i-1]*time[i-1])
		}
		return(list(auc=auc, aumc=aumc))
	}
	
	# remove missing values
	data <- na.omit(data.frame(conc, time))
	data <- data[order(data$time),]
	time <- data$time
	conc <- data$conc
	n <- nrow(data) 	           

	# subtraction of pre dosing concentration for single dose studies
	if (prev > 0) {conc <- conc - prev}

	# remove values below zero
	if (any(conc < 0)) {
		for (i in 1:length(conc)) {
			if (conc[i] < 0) {conc[i] <- NA}
		}
		warning('concentration below zero were omitted')
		data <- na.omit(data.frame(conc, time))
		time <- data$time
		conc <- data$conc	
		n <- nrow(data) 
	}

	# check input parameters
	if (numtail < 2) {stop('number of points for tail area correction must be greater than 1')}
	if (numintp < 2) {stop('number of points for interpolation must be greater than 1')}
	         
	# calculate observed auc and aumc
	auc.observed  <- add(time=time, conc=conc)$auc        
	aumc.observed <- add(time=time, conc=conc)$aumc
    
	# calculate auc from 0 to infinity and aumc from 0 to infinity  by using last numtail points above zero
	tail <- subset(data.frame(conc, time), conc > 0)
	tail <- tail[(nrow(tail)-numtail+1) : nrow(tail), ]	
	lamda <- as.real(lm(log(tail$conc)~tail$time)$coef[2])*(-1)
	auc.infinity <- auc.observed + conc[n]/lamda 
	aumc.infinity <- aumc.observed + (conc[n]*time[n])/lamda + conc[n]/lamda**2
	if(lamda < 0){
		warning('tail area correction incorrect due to increasing concentration of last numtail points')	
		auc.infinity <- NA
		aumc.infinity <- NA
	}

	# calculate auc and aumc from 0 to exact where exact must be greater than time[n-1]
	auc.interpol <- NA; aumc.interpol <- NA	
	if (!is.na(exact) & exact > time[n-1] & exact == time[n]) { # special case
		auc.interpol <- auc.observed; aumc.interpol <- aumc.observed
	}	
	if (!is.na(exact) & exact > time[n-1] & exact != time[n]) {
		conc[n] <- linpol(conc=conc[(n-numintp+1):n], time=time[(n-numintp+1):n], exact=exact)
		time[n] <- exact
		if(conc[n] < 0){warning('interpolated value below zero')}
		auc.interpol <- add(time=time, conc=conc)$auc
		aumc.interpol <- add(time=time, conc=conc)$aumc		
	} 
	

	# define output object
	res <- data.frame(AUC=c(auc.observed, auc.interpol, auc.infinity), 
		AUMC=c(aumc.observed, aumc.interpol, aumc.infinity))
	rownames(res) <- c('observed', 'interpolated', 'infinity')
	return(res)      
}      
  


# biexponential model to estimate inital and terminal half-life
biexp <- function(conc, time, prev=0, tol=1E-9){

	# get start values for optim by curve peeling
	curve.peeling <- function(x, y){

		n <- length(y)	
		res <- NA
		Fmin <- Inf
		
		for(i in 3:(n-3)) {
			parms <- lm(log(y[(i+1):n])~x[(i+1):n])
			b2 <- parms$coef[2]
			a2 <- exp(parms$coef[1])
			ynew <- abs(y - a2*exp(b2*x))
			if(!any(is.na(ynew)) && !any(ynew==0) && !any(ynew==Inf)){
				parms <- lm(log(ynew[1:i])~x[1:i])
				b1 <- parms$coef[2]
				a1 <- exp(parms$coef[1])
				F <- sum((y-(a1*exp(b1*x)+a2*exp(b2*x)))*(y-(a1*exp(b1*x)+a2*exp(b2*x))))
				b1 <- b1*(-1)
				b2 <- b2*(-1)
				if (!is.na(F) && F < Fmin && !is.na(all(b1,b2)) && all(b1>0,b2>0) && b1 > b2) {
					res <- as.real(c(a1=a1, dl=log(b1)-log(b2), a2=a2, b2=log(b2)))
					Fmin <- F
				}	
			}
		}

		if (is.na(any(res))){
			parms <- lm(log(y)~x)
			b <- parms$coef[2]
			a <- exp(parms$coef[1])
			b <- b*(-1)
			F <- sum(parms$resid*parms$resid)
			if (!is.na(F) && F < Fmin && !is.na(b) && b > 0){
				res <- as.real(c(a=a, b=log(b)))
				Fmin <- F
			}
		}

		return(res)
	}

	# remove missing values
	data <- na.omit(data.frame(conc, time))
	time <- data$time
	conc <- data$conc		

	# subtraction of pre administration concentration for single dose studies
	if (prev > 0) {conc <- conc - prev}

	# remove values below or equal to zero
	if (any(conc <= 0)) {
		for (i in 1:length(conc)) {
			if (conc[i] <= 0) {conc[i] <- NA}
		}
		warning('concentration below or equal to zero were omitted')
		data <- na.omit(data.frame(conc, time))
		conc <- data$conc
		time <- data$time		
	}
	if (nrow(data) < 4) {stop('a minimum of 4 observations are required')}

	biexploss <- function(par){
		a1 <- par[1]
		dl <- par[2]
		a2 <- par[3]
		b2 <- par[4]
		sum((conc - a1*exp(-(exp(b2) + exp(dl))*time) - a2*exp(-exp(b2)*time))^2)
	}

	singleloss <- function(par){
		a <- par[1]
		b <- par[2]
		sum((conc - a*exp(-exp(b)*time))^2)
	}

	start <- curve.peeling(y=conc, x=time)
	type <- as.character(length(start))	

	if(type == "4" && sum((conc - start[1]*exp(-(exp(start[4]) + exp(start[2]))*time) - 
			start[3]*exp(-exp(start[4])*time))^2 == Inf)) {type <- "1"}

	if(type == "2" && sum((conc - start[1]*exp(-exp(start[2])*time))^2 == Inf)){type <- "1"}

	switch(type, 
		"4" = {sol <- optim(par=start, fn=biexploss,  method=c("Nelder-Mead"), control=list(reltol=tol))	
			b1 <- (exp(sol$par[4]) + exp(sol$par[2]))
			a1 <- sol$par[1]
			b2 <- exp(sol$par[4])
			a2 <- sol$par[3]},
		"2" = {sol <- optim(par=start, fn=singleloss,  method=c("Nelder-Mead"), control=list(reltol=tol))
			b1 <- exp(sol$par[2]) 
			a1 <- sol$par[1]
			b2 <- exp(sol$par[2])
			a2 <- sol$par[1]},
		"1" = {a1 <- NA; b1 <- NA; a2<- NA; b2 <- NA},
	)

	# calculate halflife
	init.hl <- log(2) / b1
	term.hl <- log(2) / b2

	# format output object
	parms <- data.frame(initial=as.real(c(init.hl, b1, a1)),
			terminal=as.real(c(term.hl, b2, a2)))
	rownames(parms) <- c('halflife', 'slope', 'intercept')
	res <- list(parms=parms, conc=conc, time=time, method="biexp")
	class(res) <- 'halflife'
	return(res)
}

# function for halflife estimation according to the method of Lee et al.
lee <- function(time, conc, points=3, prev=0, method=c("lad", "ols", "hub", "npr")) {

	# function for lad regression
	lad <- function(y, x) {
		resid <- Inf
		for (i in 1:length(y)) {
			for (j in 1:length(y)) {
				if ((i != j) & (x[j] != x[i])) {
					slope <- (y[j] - y[i]) / (x[j] - x[i])
					intct <- y[j] - slope*x[j] 
					absresid <- abs(y - (intct + slope*x))
					if (sum(absresid) < resid) {
						mad <- median(absresid)
						resid <- sum(absresid)
						k <- slope
						d <- intct
					}
				}
			}
		}

		return(list(k=k, d=d, resid=resid, mad=mad))
	}


	# function for huber m regression 
	# acknowledgment to werner engl
	hub <- function(y, x, mad, sigmafactor=1.483, kfactor=1.5) { 

		hubloss <- function(kd) { # Huber loss for k=kd[1], d=kd[2]
			absresid <- abs(y-kd[[1]]*x-kd[[2]])
			khuber <- kfactor*sigmafactor*mad
			sum(ifelse(absresid < khuber, absresid*absresid, khuber*(2*absresid-khuber))) 
		}

		start <- as.vector(c(lm(y~x)$coef[2], lm(y~x)$coef[1]))
		res <- optim(start, hubloss, method="Nelder-Mead", control=c(reltol=1e-9))		
		return(list(k = res$par[1], d = res$par[2], resid = res$value))		
	}

	# function for nonparametric regression
	npr <- function(y, x) { 

		weighted.median <- function(w, x) { 
			data <- data.frame(x, w)
			data <- data[order(data$x),]
			i <- 1; while(sum(data$w[1:i]) < 0.5) {i <- i + 1}
			ifelse (sum(data$w[1:i-1]) == 0.5, return((data$x[i-1]+data$x[i])/2), return(data$x[i]))
		}

		total <- 0
		for(i in 1:(length(y)-1)) {
			for(j in (i+1):length(y)){total <- total + abs(x[i]-x[j])}
		}

		l <- 1
		b <- array(1:(length(y)*(length(y)-1)/2))
		w <- array(1:(length(y)*(length(y)-1)/2))
		for(i in 1:(length(y)-1)) {
			for(j in (i+1):length(y)){	
				b[l] <- ((y[i]-y[j])/(x[i]-x[j]))
				w[l] <- abs(x[i]-x[j]) / total			 
				l <- l + 1									
			}
		}

		data <- subset(data.frame(w=as.vector(w), b=as.vector(b)), b != Inf & b != -Inf)
		k <- weighted.median(w=data$w, x=data$b)
		d <- median(y-k*x)
		e <- y-k*x-d
		resid <- sum((rank(e) - 1/2*(length(y)+1))*e)
		return(list(k=k, d=d, resid=resid))
	}	

	# function for internal ols regression
	ols <- function(y, x){
		res <- lm(y~x)
		return(list(k = res$coef[2], 
		    	d = res$coef[1], 
		    	resid = sum(res$resid*res$resid)))
	}

	# exclude missing values
	data <- na.omit(data.frame(conc, time))

	# subtraction of pre administration concentration for single dose studies
	if (prev > 0) {data$conc <- data$conc - prev}

	# check input parameters and remove values below or equal to zero
        method = match.arg(method)
	if (any(data$time < 0)) {stop('timepoint below zero')}
	if (points < 2) {stop('not enough points in terminal phase')}

	# remove values below or equal to zero
	if (any(data$conc <= 0)) {
		for (i in 1:nrow(data)) {
			if (data$conc[i] <= 0) {data$conc[i] <- NA}
		}
		warning('concentration below or equal to zero were omitted')
		data <- na.omit(data)	
	}
	if (nrow(data) < 4) {stop('a minimum of 4 observations are required')}
  
	# transform data by logarithm at base 10
	n <- nrow(data)
 	data$conc <- log10(data$conc)
	conc <- data$conc
	time <- data$time
	
	# calculate parameters of one-phase model
	switch(method, 
		"lad"={model <- lad(y=conc, x=time)}, 
		"ols"={model <- ols(y=conc, x=time)}, 
		"hub"={model <- hub(y=conc, x=time, mad=lad(y=conc, x=time)$mad)},
		"npr"={model <- npr(y=conc, x=time)
	},)

	# inital halflife = terminal halflife for one-phase model
	final.term.model <- model
	final.init.model <- model
	resid <- model$resid
	final.chgpt <- NA

	# check special cases
	if(model$k >= 0) {
		resid <- Inf
		final.term.model$k <- NA
		final.init.model$k <- NA
	}
	
	
	# calculate parameters of two-phase models
	if (points > n-2) {stop('not enough points for inital phase')}
	for (i in 2:(n-points)) {

		init.conc <- conc[1:i]
		init.time <- time[1:i]
		term.conc <- conc[(i+1):n]
		term.time <- time[(i+1):n]

		# calculate parameters of two-phase model
		switch(method, 
			"lad"={
				init.model <- lad(y=init.conc, x=init.time)
				term.model <- lad(y=term.conc, x=term.time)
		},	"ols"={
				init.model <- ols(y=init.conc, x=init.time)
				term.model <- ols(y=term.conc, x=term.time)
		}, 	"hub"={
				init.model <- hub(y=init.conc, x=init.time, mad=lad(y=init.conc, x=init.time)$mad)
				term.model <- hub(y=term.conc, x=term.time, mad=lad(y=term.conc, x=term.time)$mad)
		}, 	"npr"={
				init.model <- npr(y=init.conc, x=init.time)
				term.model <- npr(y=term.conc, x=term.time)
		},)


		# check changeover criteria and for negative slopes 
		lower <- data$time[i]
		upper <- data$time[i+1]
		chgpt <- (init.model$d - term.model$d) / (term.model$k - init.model$k)
		if (!(chgpt <= lower | chgpt >= upper) & 
			(term.model$k < 0) & 
			(init.model$k < 0) & 
			(init.model$k <= term.model$k)){  
  			if (sum(term.model$resid, init.model$resid) < resid) {
				final.init.model <- init.model
				final.term.model <- term.model
				final.chgpt <- as.real(chgpt)
				resid <- sum(term.model$resid, init.model$resid)
			}  
		}
	}	

	init.hl <- -log10(2)/final.init.model$k
	term.hl <- -log10(2)/final.term.model$k

	# format output objects
	parms <- data.frame(initial=as.real(c(init.hl, final.init.model$k, final.init.model$d)),
			terminal=as.real(c(term.hl, final.term.model$k, final.term.model$d)))
	rownames(parms) <- c('halflife', 'slope', 'intercept')
	res <- list(parms=parms, chgpt=as.real(final.chgpt), conc=10**conc, time=time, method='lee')
	class(res) <- 'halflife'
	return(res)
	
}



# plot function for halflife objects 
plot.halflife <- function(x, xlab='Time', ylab='Concentration', main='Half-life Estimation', xlim=NULL, ylim=NULL, ...) {

	init.k <- x$parms[2,1] 
	init.d <- x$parms[3,1] 
	term.k <- x$parms[2,2] 
	term.d <- x$parms[3,2] 

	if(is.null(xlim)){xlim <- c(min(x$time), max(x$time))}
	if(is.null(ylim)){ylim <- c(min(x$conc), max(x$conc))}

	plot(x$time, x$conc, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, main=main, ...)

	switch(x$method,
		"lee"={
			if (is.na(x$chgpt)) {x$chgpt <- min(x$time)}
			plot(function(x) 10**(init.k*x+init.d), xlim[1], x$chgpt, add=TRUE)
			plot(function(x) 10**(term.k*x+term.d), x$chgpt, xlim[2], add=TRUE) 	
	},	"biexp"={
			if (init.k == term.k){
				plot(function(x) init.d*exp(-init.k*x), xlim[1], xlim[2], add=TRUE)
			}
			if (init.k != term.k){
				plot(function(x) init.d*exp(-init.k*x)+term.d*exp(-term.k*x), xlim[1], xlim[2], add=TRUE)
			}
	},)

		
}	

