###### define lin pred. of various ######
###### dose-reponse models         ######
#########################################
pow = function(dose,p=1,dmax=NULL,dc=NULL,off=NULL,scale=NULL) {
    if (is.null(off)) off=(p<=0)
    if (is.null(scale)) scale=1

    if (length(p)==1) {
        dose=dose/scale+off
        if (p==0) {
            d=cbind(log(dose));colnames(d)="log.dose"} else {    
            d=cbind(dose^p); colnames(d)=paste("dose",ifelse(p==1,"",paste("^",p,sep="")), sep="")
        }
    attr(d,"power")=p; attr(d,"off")=off; attr(d,"scale")=scale
    return(d)
    }
    
    if (is.null(dmax)) dmax=max(dose)/scale+off else dmax=dmax/scale+off  #default for dmax is maximum dose
    if (length(p)==2 && all(p!=0) ) {
        d.pow1=(dose/scale+off[1])^p[1]; d.pow2=(dose/scale+off[2])^p[2]
        d=cbind(d.pow1, d.pow2);
        colnames(d)<- c(paste("dose",ifelse(p[1]==1,"",paste("^",p[1],sep="")), sep=""),
                        paste("dose",ifelse(p[2]==1,"",paste("^",p[2],sep="")), sep=""))
        attr(d,"power")=p; attr(d,"dmax")=dmax; attr(d,"off")=off; attr(d,"scale")=scale
        return(d)
    }

    if (length(p)==2 && any(p==0) ) {
        off=c(off[match(0,p)],off[-match(0,p)]); dmax=c(dmax[match(0,p)],dmax[-match(0,p)])
        p=c(0,p[-match(0,p)]);  # log is always at power[1]
        log.d=log(dose/scale+off[1]); d.pow=(dose/scale+off[2])^p[2]
        d=cbind(log.d, d.pow);
        colnames(d)<- c("log.dose",paste("dose",ifelse(p[2]==1,"",paste("^",p[2],sep="")), sep=""))
        attr(d,"power")=p; attr(d,"dmax")=dmax; attr(d,"off")=off; attr(d,"scale")=scale
        return(d)
    }       
    
    if (length(p)==3 && all(p!=0) ) {
        d.pow1=(dose/scale+off[1])^p[1]; d.pow2=(dose/scale+off[2])^p[2]; d.pow3=(dose/scale+off[3])^p[3]
        d=cbind(d.pow1, d.pow2, d.pow3)
        colnames(d)<- c(paste("dose",ifelse(p[1]==1,"",paste("^",p[1],sep="")), sep=""),
                        paste("dose",ifelse(p[2]==1,"",paste("^",p[2],sep="")), sep=""),
                        paste("dose",ifelse(p[3]==1,"",paste("^",p[3],sep="")), sep=""))
        attr(d,"power")=p; attr(d,"dmax")=dmax; attr(d,"dc")=dc+off; attr(d,"off")=off; attr(d,"scale")=scale
        return(d)
    }

    if (length(p)==3 && p[1]==0 ) {  # log is always at power[1]
        log.d=log(dose/scale+off[1]); d.pow1=(dose/scale+off[2])^p[2]; d.pow2=(dose/scale+off[3])^p[3]
        d=cbind(log.d, d.pow1, d.pow2);
        colnames(d)<- c("log.dose",paste("dose",ifelse(p[2]==1,"",paste("^",p[2],sep="")), sep=""),
                        paste("dose",ifelse(p[3]==1,"",paste("^",p[3],sep="")), sep=""))
        attr(d,"power")=p; attr(d,"dmax")=dmax; attr(d,"dc")=dc+off; attr(d,"off")=off; attr(d,"scale")=scale
        return(d)
    }


}

expo = function(dose,order=1,scale=1,off=0) {
    if(all(order!=c(1,2))) stop(paste("only exponential (order=1) and double exponential models (order=2) are supported"))
    if (order==1) {
            d=cbind(exp(dose/scale+off));colnames(d)="exp.dose"} else {  
            d=cbind(exp(exp(dose/scale+off))); colnames(d)="exp.exp.dose"
        }
    attr(d,"expo")=order; attr(d,"scale")=scale; attr(d,"off")=off;
    return(d)
    }

compart = function(dose, dmax=NULL, logd=FALSE, off=ifelse(logd,1,0), scale=1) {
    if (is.null(dmax)) dmax=max(d)
    if (logd) {d=log(dose/scale+off); dmax=log(dmax/scale+off)} else {d=dose/scale+off; dmax=dmax/scale+off}
    attr(d,"compart")=TRUE; attr(d,"dmax")=dmax; attr(d,"logd")=logd; attr(d,"scale")=scale; attr(d,"off")=off;
    return(d)
    }

emax = function(dose, ED50=NULL, logd=FALSE, off=ifelse(logd,1,0), scale=1) {
    if (is.null(ED50)) ED50=((max(dose)-min(dose))/2)
    if (logd) {d=log(dose/scale+off); ED50=log(ED50/scale+off)} else {d=dose/scale+off; ED50=ED50/scale+off} 
    attr(d,"emax")=TRUE; attr(d,"ED50")=ED50; attr(d,"logd")=logd; attr(d,"scale")=scale; attr(d,"off")=off;
    return(d)
    }

semax = function(dose, ED50=NULL, hill=NULL, logd=FALSE, off=ifelse(logd,1,0), scale=1) {
    if (is.null(ED50)) ED50=((max(dose)-min(dose))/2)
    if (logd) {d=log(dose/scale+off); ED50=log(ED50/scale+off)} else {d=dose/scale+off; ED50=ED50/scale+off} 
    if (is.null(hill)) hill=1
    attr(d,"semax")=TRUE; attr(d,"ED50")=ED50; attr(d,"hill")=hill; attr(d,"logd")=logd; attr(d,"scale")=scale; attr(d,"off")=off;
    return(d)
    }

logistic = function(dose, ED50=NULL, slope=NULL, logd=FALSE, off=ifelse(logd,1,0), scale=1) {
    if (is.null(ED50)) ED50=((max(dose)-min(dose))/2)
    if (logd) {d=log(dose/scale+off); ED50=log(ED50/scale+off)} else {d=dose/scale+off; ED50=ED50/scale+off} 
    if (is.null(slope)) slope=1
    attr(d,"logistic")=TRUE; attr(d,"ED50")=ED50; attr(d,"slope")=slope; attr(d,"logd")=logd; attr(d,"scale")=scale; attr(d,"off")=off;
    return(d)
    }

weibull = function(dose, infl=NULL, slope=NULL, logd=FALSE, off=ifelse(logd,1,0), scale=1) {
    if (is.null(infl)) infl=((max(dose)-min(dose))/2)
    if (logd) {d=log(dose/scale+off); infl=log(infl/scale+off)} else {d=dose/scale+off; infl=infl/scale+off}
    if (is.null(slope)) slope=1
    attr(d,"weibull")=TRUE; attr(d,"infl")=infl; attr(d,"slope")=slope; attr(d,"logd")=logd; attr(d,"scale")=scale; attr(d,"off")=off;
    return(d)
    }


############################################
# find initial parameter estimates ########
# for plotting candidate models     ########
############################################

getpi = function(M,d,low,high) {         
     link.low=M$family$linkfun(low)
     link.high=M$family$linkfun(high)

     if (is.null(M$model)) {M$model=d; attr(M$model,"power")=1; attr(M$model,"off")=0; attr(M$model,"scale")=1}  #if no model specified, assume linpred=dose
     ### Power Models ###
     if (!is.null(attr(M$model,"power"))) {
        power=attr(M$model,"power")
        dose=do.call("pow",list(dose=d,p=power,dmax=attr(M$model,"dmax"),
                                off=attr(M$model,"off"), scale=attr(M$model,"scale")))

        ### 1-parm ###
        if (length(power)==1) {
            beta=(link.high - link.low)/(dose[length(dose)]-dose[1])
            alpha=link.low-beta*dose[1]
            return(M$family$linkinv(alpha+beta*dose))
        }   
        ### 2-parms  ###
        dmax=attr(M$model,"dmax")
        if (length(power)==2 && all(power!=0) ) {
            A=matrix(NA,3,3); b=matrix(NA,3,1)
            A[1,]=c(1,dose[1,1],dose[1,2])
            A[2,]=c(1,dmax[1]^power[1],dmax[2]^power[2])
            A[3,]=c(0,power[1]*dmax[1]^(power[1]-1),power[2]*dmax[2]^(power[2]-1))
            b[1]=link.low
            b[2]=link.high
            b[3]=0
            b=solve(A)%*%b
            return(M$family$linkinv(b[1]+b[2]*dose[,1]+b[3]*dose[,2]))
        }

        if (length(power)==2 && any(power==0) ) {
            A=matrix(NA,3,3); b=matrix(NA,3,1)
            A[1,]=c(1,dose[1,1],dose[1,2])
            A[2,]=c(1,log(dmax[1]),dmax[2]^power[2])
            A[3,]=c(0,1/dmax[1],power[2]*dmax[2]^(power[2]-1))
            b[1]=link.low
            b[2]=link.high
            b[3]=0
            b=solve(A)%*%b
            return(M$family$linkinv(b[1]+b[2]*dose[,1]+b[3]*dose[,2]))
        }
        ### 3-parms  ###
        dmax=attr(M$model,"dmax")
        dc=attr(M$model,"dc")
        if (length(power)==3 && all(power!=0) ) {
            A=matrix(NA,4,4); b=matrix(NA,4,1)
            A[1,]=c(1,dose[1,1],dose[1,2],dose[1,3])
            A[2,]=c(1,dmax[1]^power[1],dmax[2]^power[2],dmax[3]^power[3])
            A[3,]=c(0,power[1]*dmax[1]^(power[1]-1),power[2]*dmax[2]^(power[2]-1),power[3]*dmax[3]^(power[3]-1))
            A[4,]=c(0,power[1]*(power[1]-1)*dc[1]^(power[1]-2),power[2]*(power[2]-1)*dc[2]^(power[2]-2),power[3]*(power[3]-1)*dc[3]^(power[3]-2))
            b[1]=link.low
            b[2]=link.high
            b[3]=0
            b[4]=0
            beta=solve(A)%*%b
            return(M$family$linkinv(beta[1]+beta[2]*dose[,1]+beta[3]*dose[,2]+beta[4]*dose[,3]))
        }

        if (length(power)==3 && any(power==0) ) {
            A=matrix(NA,4,4); b=matrix(NA,4,1)
            A[1,]=c(1,dose[1,1],dose[1,2],dose[1,3])
            A[2,]=c(1,log(dmax[1]),dmax[2]^power[2],dmax[3]^power[3])
            A[3,]=c(0,1/dmax[1],power[2]*dmax[2]^(power[2]-1),power[3]*dmax[3]^(power[3]-1))
            A[4,]=c(0,-1/dc[1]^2,power[2]*(power[2]-1)*dc[2]^(power[2]-2),power[3]*(power[3]-1)*dc[3]^(power[3]-2))
            b[1]=link.low
            b[2]=link.high
            b[3]=0
            b[4]=0
            beta=solve(A)%*%b
            return(M$family$linkinv(beta[1]+beta[2]*dose[,1]+beta[3]*dose[,2]+beta[4]*dose[,3]))
        }

    } #end power

    ### exp model ###
    if (!is.null(attr(M$model,"expo"))) {
        dose=do.call("expo",list(dose=d,order=attr(M$model,"expo"),
                                scale=attr(M$model,"scale"), off=attr(M$model,"off")))
        beta=(link.high - link.low)/(dose[length(dose)]-dose[1])
        alpha=link.low-beta*dose[1]
        return(M$family$linkinv(alpha+beta*dose))
    }

    ### compartment model ###  can use all possible link for plotting, but currently implemented only with logit link
    if (!is.null(attr(M$model,"compart"))) {
        dose=do.call("compart",list(dose=d, dmax=attr(M$model,"dmax"), logd=attr(M$model,"logd"),
                                scale=attr(M$model,"scale"), off=attr(M$model,"off")))
        dmax=attr(M$model,"dmax")
        alpha=link.low
        gamma=dmax
        beta=(link.high-link.low)/(dmax*exp(-1))
        return(M$family$linkinv(alpha+beta*dose*exp(-dose/gamma)))
    }

    ### emax model ###  can use all possible links for plotting, but currently implemented only with logit link
    if (!is.null(attr(M$model,"emax"))) {
        dose=do.call("emax",list(dose=d, ED50=attr(M$model,"ED50"), logd=attr(M$model,"logd"), 
                                scale=attr(M$model,"scale"), off=attr(M$model,"off")))
        ED50=attr(M$model,"ED50")
        alpha=link.low
        beta=(link.high-link.low)*(ED50+max(dose))/max(dose)
        #gamma=ED50*(2*beta/(link.high-link.low)-1)
        return(M$family$linkinv(alpha+beta*dose/(ED50+dose)))
    }

    ### semax model ###  can use all possible links for plotting, but currently implemented only with logit link
    if (!is.null(attr(M$model,"semax"))) {
        dose=do.call("semax",list(dose=d, ED50=attr(M$model,"ED50"), hill=attr(M$model,"hill"), logd=attr(M$model,"logd"),
                                scale=attr(M$model,"scale"), off=attr(M$model,"off")))
        ED50=attr(M$model,"ED50")
        hill=attr(M$model,"hill")
        alpha=link.low
        beta=(link.high-link.low)*(ED50^hill+max(dose)^hill)/max(dose)^hill
        #beta=(max(dose)^hill-ED50^hill)*(link.high-link.low)/(max(dose)^hill-2*ED50^hill)
        #gamma=(ED50^hill*(2*beta/(link.high-link.low)-1))^(1/hill)
        return(M$family$linkinv(alpha+beta*dose^hill/(ED50^hill+dose^hill)))
    }

    ### 4-parameter logistic model ###  can use all possible links for plotting, but currently implemented only with logit link
    if (!is.null(attr(M$model,"logistic"))) {
        dose=do.call("logistic",list(dose=d, ED50=attr(M$model,"ED50"), slope=attr(M$model,"slope"), logd=attr(M$model,"logd"),
                                scale=attr(M$model,"scale"), off=attr(M$model,"off")))
        ED50=attr(M$model,"ED50")
        slope=attr(M$model,"slope")
        beta=(link.high-link.low)*(1/(1+exp((ED50-max(dose))/slope)) + 1/(1+exp(ED50/slope))) 
        alpha=link.low-beta/(1+exp(ED50/slope))
        #beta=(max(dose)^hill-ED50^hill)*(link.high-link.low)/(max(dose)^hill-2*ED50^hill)
        #gamma=(ED50^hill*(2*beta/(link.high-link.low)-1))^(1/hill)
        return(M$family$linkinv(alpha+beta/(1+exp((ED50-dose)/slope))))
        }

    ### 4-parameter weibull model ###  can use all possible links for plotting, but currently implemented only with logit link
    if (!is.null(attr(M$model,"weibull"))) {
        dose=do.call("weibull",list(dose=d, infl=attr(M$model,"infl"), slope=attr(M$model,"slope"), logd=attr(M$model,"logd"),
                                scale=attr(M$model,"scale"), off=attr(M$model,"off")))
        gamma=attr(M$model,"infl")
        slope=attr(M$model,"slope")
        beta=(link.high-link.low)/(exp((-exp(slope*(max(dose)-gamma)))) - exp((-exp(-slope*gamma)))) 
        alpha=link.low-beta*exp((-exp(-slope*gamma)))
        #beta=(max(dose)^hill-ED50^hill)*(link.high-link.low)/(max(dose)^hill-2*ED50^hill)
        #gamma=(ED50^hill*(2*beta/(link.high-link.low)-1))^(1/hill)
        return(M$family$linkinv(alpha + beta*exp((-exp(slope*(dose-gamma)))) ))
        }

}


############################################
#### Plots candidate models             ####
#### for given placebo effect (low)     ####
#### and maximum efficacy (high)        ####
############################################

plotModels=function(dose,models,low=0.0001,high=0.9999,label=names(models), steps=min(dose[-1])/100, print=FALSE, ...)
{
low=ifelse(low==0,0.0001,low);high=ifelse(high==1,0.9999,high);
if (low>=high) stop(paste("0<low<high<1"))
m=length(models)   #number of cand. models 
d=seq(min(dose),max(dose),by=steps) 
d=sort(unique(c(d,dose)))
pi=matrix(-99,length(d),m+1)
pi[,1]=d

for (i in 1:m)
    {
    pi[,i+1]=getpi(models[[i]],d,low,high)
    }

########################################
########   Plotting Part ###############
########################################
library(lattice)
if (is.null(label)) label=paste("M", 1:m, sep = "")
model.names=factor(1:m,labels=label)
plotData = data.frame(doses = rep(pi[,1],m), resp = matrix(pi[,-1],ncol=1), model = rep(model.names, each=length(pi[,1]))) 
if (print) return(plotData)    # if print=TRUE return the matrix used to plot data without the plot, otherwise, plot

plot.symbol = trellis.par.get("plot.symbol")
plot.symbol$pch = 15
plot.symbol$cex=0.8
trellis.par.set("plot.symbol", plot.symbol)

xyplot(resp~doses|model, data = plotData, as.table = TRUE,
                        xlab="dose", ylab="efficacy",  
            panelData = data.frame(d = dose),
            panel = function(x, y,...,panelData)
                {
                panel.abline(h = c(low, high), lty = 2)
                panel.xyplot(x, y, type="l")
                ind <- match(panelData$d, x)
                panel.xyplot(x[ind], y[ind]) # plot squares
                }
       ,...)

}

### Functions for defining models
#emax  not imlemented yet

#myLatticeSettings <- function()
#    list(layout.heights =
#         list(top.padding = 0,
#              main.key.padding = 0,
#              key.axis.padding = 0,
#              axis.xlab.padding = 0,
#              xlab.key.padding = 0,
#              key.sub.padding = 0,
#              bottom.padding = 0),
#         layout.widths =
#         list(left.padding = 0,
#              key.ylab.padding = 0,
#              ylab.axis.padding = 0,
#              axis.key.padding = 0,
#              right.padding = 0)
#         )
#
#trellis.par.set(myLatticeSettings())