R/plot.R

In DavidQuigley/HTDoseResponseCurve: Dose Response Curve Evaluation from Incucyte and Other High-Throughput Platforms

#' Plot a grid of colors given a matrix of numbers
#' 
#' @param M matrix of values to plot
#' @param block.height element height, default 20
#' @param block.width element width, default 10
#' @param space.X space between elements on horizontal axis, default 3
#' @param space.Y space between elements on vertical axis, default 10
#' @param cex.x cex value for x axis
#' @param cex.y cex value for y axis
#' @param border boolean, draw a border around grid elements
#' @param color_palatte vector of colors c(cold,middle,hot), defaults to 
#' c("blue", "white", "red") 
#' @param color_bounds two element vector with lower and upper limit for color 
#' palatte, defaults to c( min(M), max(M) )
#' @param num_colors number of color gradations in palatte, default 50
#' @param main Plot main title, default blank.
#' @param xlab Plot x axis label, default blank.
#' @param ylab Plot y axis label, default blank.
#' @param plot_values draw element values in the grid blocks, default FALSE
#' @return M, including possible truncation from color bounds
#' @examples 
#' M = matrix(1:12, nrow=3, ncol=4,byrow=TRUE)
#' plot_color_grid(M)
#' dimnames(M)[[1]] = c("row_a","row_b","row_c")
#' dimnames(M)[[2]] = c("col1", "col2", "col3","col4")
#' plot_color_grid(M, color_palatte=c("blue","white","yellow" ) )
#' @importFrom graphics plot rect axis lines points text
#' @importFrom grDevices colorRampPalette
#' @export
plot_color_grid = function(M, block.height=20, block.width=10, space.X=3, 
                         space.Y=10, cex.x=1, cex.y=1, border=TRUE, 
                         color_palatte=c("blue","white","red"), 
                         color_bounds=NA, num_colors=50, 
                         main="", xlab="", ylab="", plot_values=FALSE ) {
    
    cmap = grDevices::colorRampPalette( colors=color_palatte )( num_colors )
    if( is.na(color_bounds)[1] ){
        color_bounds = c( min(M, na.rm=TRUE), max(M, na.rm=TRUE) )
    }
    if( sum(M>color_bounds[2], na.rm=TRUE)>0 ){
        warning("value(s) in M exceed maximum color bound, truncating value(s)")
    }
    if( sum(M<color_bounds[1], na.rm=TRUE)>0 ){
        warning("value(s) in M exceed minimum color bound, truncating value(s)")
    }
    M[M>color_bounds[2]] = color_bounds[2]
    M[M<color_bounds[1]] = color_bounds[1]
    MC = color_scale( M, cmap, color_bounds = color_bounds )
    
    n.rows = dim(MC)[1]
    n.cols = dim(MC)[2]
    total.width =  ( block.width*n.cols) + ( (n.cols-1) * space.X)
    total.height = ( block.height * n.rows ) + ( (n.rows-1) * space.Y)
    graphics::plot(0, 0, col="white", xlim=c(0,total.width), 
                   ylim=c(0,total.height), axes=FALSE, xlab=xlab, ylab=ylab, 
                   bg="azure2", main=main)
    xlab_locs = rep(0, n.cols)
    ylab_locs = rep(0, n.rows)
    cur.y = total.height
    for(rr in 1:n.rows){
        for(cc in 1:n.cols){  
            this.x.left = (cc-1)*block.width + (cc-1)*space.X
            this.x.right = this.x.left+block.width
            xlab_locs[cc] = this.x.right - (block.width)
            if( border ){
                graphics::rect( this.x.left , cur.y - block.height, 
                            this.x.right, cur.y, col=MC[rr,cc], border="azure2")
            }else{
                graphics::rect( this.x.left , cur.y - block.height,this.x.right, 
                                cur.y, col=MC[rr,cc], border=NA)
            }
            if( plot_values ){
                graphics::text( this.x.right - (block.width/2), 
                      cur.y - (block.height/2), 
                      M[rr,cc], cex=0.75 )
            }
        }
        ylab_locs[rr] = cur.y - (0.5*block.height)
        cur.y = cur.y - block.height - space.Y
    }
    graphics::axis(1, at=xlab_locs, labels=dimnames(MC)[[2]], las=2, 
                   cex.axis=cex.x, tick=FALSE, padj=1, line=-1 )
    graphics::axis(2, at=ylab_locs, labels=dimnames(MC)[[1]], las=2, 
                   cex.axis=cex.y, tick=FALSE, hadj=1, line=-1.5)
    M
}

#' Report colors that will be used when plotting a HTfit object. 
#' 
#' If a vector of colors is passed to col, these colors will be used instead 
#' of the default colors. 
#' 
#' @param F HTfit object
#' @param col optional vector of valid R colors to be used for the plot; 
#' defaults to the color scheme built into HTDoseResponseCurve.
#' @return data frame of unique condition, color
#' @export
HT_fit_plot_colors=function( F, col=NULL ){
    N = length(F$unique_conditions)
    if( is.null(col) ){
        col = INC_colors_DRC[1:N]
    }else{
         if(length(col)<N){
            stop(paste("Must provide",N,"colors to plot this figure"))   
         }
    }
    uc = as.character( unique( F$input$conditions_to_fit ) )
    
    data.frame(condition=uc, col=col, stringsAsFactors=FALSE)
}

#' Plot a dose response curve from a HTfit object
#' 
#' If the curve cannot be fit (meaning that the optimization method used by 
#' the drm() method from the drc package failes to converge) then the summarized 
#' points will be plotted without a curve connecting them.
#' 
#' @param x HTfit object
#' @param ... standard parameters for \code{plot} function
#' @param bar_multiple multiplier for standard error bars, default 2
#' @param summary_method summary method for points to plot in timecourse, one 
#' of ("mean", "median"), defaults "mean"
#' @param show_EC50 show dotted vertical lines at EC50 values if they are fit,
#' default FALSE
#' @param show_SF50 show dotted vertical lines at SF50 values if they are fit,
#' default FALSE. Remember that SF50 is the dose that produces a surviving 
#' fraction of 50%, whereas EC50 is the dose that produces a drop of 50% from 
#' the minimum to maximum measured effect. Has no effect if the SF50 is not 
#' reached.
#' @param show_legend show a legend, default FALSE
#' @param show_x_log_tics if axes==TRUE and this parameter is TRUE, draw tics 
#' between each power of 10 scaled appropriately. If axes==TRUE and this 
#' parameter is FALSE, omit ticks between powers of 10. Default TRUE.
#' @param show_x_exponent If axes==TRUE and this parameter is TRUE, show labels 
#' on X axis as 10^2, 10^3, ... . If axes==TRUE and this parameter is FALSE, 
#' show labels on X axis as an exponent (2, 3, ... ). Default TRUE.
#' @param log_axis_x logical. If TRUE, X axis is to be logarithmic; if FALSE, 
#' X axis is original scale. Default is TRUE. 
#' @return none
#' @examples 
#' sample_types = rep( c(rep("line1",3), rep("line2",3)), 5)
#' treatments = c(rep("DMSO",6), rep("drug",24))
#' concentrations = c( rep(0,6),rep(200,6), rep(500,6),rep(1000,6),rep(5000,6))
#' values=c(100,99,100,90,91,92,99,97,99,89,87,88,86,89,88,56,59,58,66,65,67,
#'          25,23,24,42,43,46,4,5,9)
#' hours = rep(48, length(values))
#' plate_id = "plate_1"
#' ds = create_dataset( sample_types=sample_types, treatments=treatments, 
#'                      concentrations=concentrations, hours=hours, 
#'                      values=values, plate_id=plate_id, 
#'                      negative_control = "DMSO")
#' library(drc)
#' # Fit model using three-parameter log-logistic function
#' fit_1=fit_DRC(ds, sample_types=c("line1", "line2"), treatments=c("drug"), 
#'         hour = 48, fct=drc::LL.3() )
#' plot(fit_1)
#' plot(fit_1, show_EC50=FALSE, show_legend=FALSE, lwd=3,col=c("black", "gold"),
#'      xlim=c(0, 1e4), xlab="concentration nM", 
#'      ylab="surviving fraction")
#' legend(1.5, 0.3, c("Line 1", "Line 2"), col=c("black", "gold"), pch=15)
#' @seealso \code{\link{fit_DRC}}
#' @export
plot.HT_fit = function( x, ..., bar_multiple=2, 
                        summary_method="mean", show_EC50=FALSE, 
                        show_SF50=FALSE,
                        show_legend=FALSE, show_x_log_tics=TRUE, 
                        show_x_exponent=TRUE, log_axis_x=TRUE){
    if( summary_method!="mean" & summary_method !="median"){
        stop("parameter summary_method must be one of {mean, median}")   
    }
    
    N = length(x$unique_conditions)
    plot_parameters = list(...)
    if( length( which( names(plot_parameters)== "pch"  ) )==0 )
        plot_parameters[["pch"]] = rep(15, N)
    if( length( which( names(plot_parameters)== "col"  ) )==0 ){
        plot_parameters[["col"]] = INC_colors_DRC[1:N]
    }else{
        if(length(plot_parameters[["col"]])==1)
            plot_parameters[["col"]] = rep(plot_parameters[["col"]], N)
    }
    if( length( which( names(plot_parameters)== "lty"  ) )==0 )
        plot_parameters[["lty"]] = rep(1, N)
    if( length( which( names(plot_parameters)== "ylim"  ) )==0 )
        plot_parameters[["ylim"]] = c(0, 1.2)
    if( length( which( names(plot_parameters)== "xlab"  ) )==0 )
        plot_parameters[["xlab"]] = ""
    if( length( which( names(plot_parameters)== "ylab"  ) )==0 )
        plot_parameters[["ylab"]] = ""
    if( length( which( names(plot_parameters)== "cex"  ) )==0 )
        plot_parameters[["cex"]] = 1
    if( length( which( names(plot_parameters) == "cex.axis"))==0 )
        plot_parameters[["cex.axis"]] = 1
    if( length( which( names(plot_parameters) == "cex.lab"))==0 )
        plot_parameters[["cex.lab"]] = 1.5
    if( length( which( names(plot_parameters) == "xlim"))==0 )
        plot_parameters[["xlim"]] = c(0, 1e4)
    plot_parameters[["yaxs"]] = "i"
    plot_parameters[["xaxs"]] = "i"
    if( !log_axis_x ){
        plot_parameters[["log"]] = ""   
    }
    pc = HT_fit_plot_colors( x, col=plot_parameters[["col"]][ 1:N ] )
    cond2color = hsh_from_vectors(pc$condition, pc$col)
    
    Mstat = plyr::ddply( x$input, c("conditions_to_fit", "concentration"), 
                         function(po){ data.frame( 
                             mu=mean(po$value, na.rm=TRUE), 
                             med=stats::median(po$value, na.rm=TRUE),
                             sterr = se(po$value), 
                             bar_width=(po$concentration/10), 
                             stringsAsFactors=FALSE ) }
    )
    if( summary_method=="mean"){
        Mstat$value = Mstat$mu
    }
    else{
        Mstat$value = Mstat$med
    }
    
    # if user wants axes (either by passing axes=TRUE or not specifying), we're 
    # going to draw our own axes, so set the plot parameter axes to FALSE. If 
    # the user passed axes=FALSE and therefore DOESN'T want axes, respect that.
    if( length( which( names(plot_parameters)== "axes"  ) )==0 ){
        plot_parameters[["axes"]] = FALSE
        draw_axes = TRUE
    }else{
        draw_axes = plot_parameters[["axes"]]
        if( plot_parameters[["axes"]] ){
            plot_parameters[["axes"]] = FALSE
        }
    }
    
    if( x$is_fitted ){
        plot_parameters[["type"]] = "none"
        plot_parameters[["legend"]] = FALSE # for plot.drc()
        plot_parameters[["bp"]] = 1 # for plot.drc()
        
        plot_parameters = c( list(x$model), plot_parameters)
        do.call( graphics::plot, plot_parameters )
        graphics::box(col="white", lwd=4)
        if( show_EC50 ){
            EC50 = x$fit_stats$coef_EC50 
            for( i in 1:length(EC50) ){
                if( !is.na(EC50[i]) & !is.infinite(EC50[i]) ){
                    y_pred=stats::predict( x$model, 
                                      data.frame( concentration=EC50[i],
                                      conditions_to_fit=x$unique_conditions[i]))
                    graphics::lines( c( EC50[i], EC50[i] ), c(0, y_pred), 
                           col=hsh_get( cond2color,
                                        rownames(x$fit_stats)[i] ),
                           lwd=plot_parameters[["lwd"]], lty=3 )
                }
            }
        }
        if( show_SF50 ){
            SF50 = x$fit_stats$SF50
            for( i in 1:length(SF50) ){
                if( !is.na(SF50[i]) & !is.infinite(SF50[i]) ){
                    graphics::lines( c( SF50[i], SF50[i] ), c(0, 0.5), 
                                     col=hsh_get( cond2color,
                                                  rownames(x$fit_stats)[i] ),
                                     lwd=plot_parameters[["lwd"]], lty=3 )
                }
            }
        }
        bandwidth_factor=10
        x_legend=10
    }else{
        Mstat$concentration = log10( Mstat$concentration )
        Mstat$concentration[is.infinite( Mstat$concentration)] = 0
        conditions = unique(Mstat$conditions_to_fit)
        unfit = unique(Mstat[,c("conditions_to_fit", "concentration","value")])
        xlim_transformed=log10( plot_parameters[["xlim"]] )
        xlim_transformed[is.infinite(xlim_transformed)] = 0
        plot( -1, -1, 
        #plot( Mstat$concentration[Mstat$conditions_to_fit==conditions[1]], 
        #      Mstat$value[Mstat$conditions_to_fit==conditions[1]],
              axes=FALSE, 
              xlim=xlim_transformed,
              ylim=plot_parameters[["ylim"]],
              pch=plot_parameters[["pch"]][1],
              col=plot_parameters[["col"]][1],
              xlab=plot_parameters[["xlab"]],
              ylab=plot_parameters[["ylab"]],
              main=plot_parameters[["main"]],
              yaxs = "i",
              xaxs= "i")
        #plot_parameters[["x"]] = -1
        #do.call( graphics::plot, plot_parameters)
        Mstat$bar_width = 0.2
        bandwidth_factor=50
        x_legend=0.2
    }
    if( draw_axes ){
        graphics::axis( 2, at=c(0, 0.5, 1), 
              labels=c("0", "0.5", "1"), las=1, 
              cex.axis=plot_parameters[["cex.axis"]],
              font=plot_parameters[["font"]])
        if( log_axis_x ){
            log10_low = log10( plot_parameters[["xlim"]][1] )
            log10_high = log10( plot_parameters[["xlim"]][2] )
            if( plot_parameters[["xlim"]][1]==0 )
                log10_low=0
            tics = logtics( log10_low, log10_high,
                            show_x_log_tics, show_x_exponent)
            if( !x$is_fitted)
                tics$values = log10(tics$values)
            graphics::axis(1, at=tics$values, labels=tics$labels, las=1, 
                 font=plot_parameters[["font"]],
                 cex.axis=plot_parameters[["cex.axis"]], lwd.ticks=1)   
            graphics::axis(1, at=tics$values[tics$lwd==2],
                 labels=rep("", sum(tics$lwd==2)), las=1, 
                 font=plot_parameters[["font"]],
                 cex.axis=plot_parameters[["cex.axis"]], lwd.ticks=2)    
        }else{
            xlow = plot_parameters[["xlim"]][1]
            xhigh = plot_parameters[["xlim"]][2]
            values = seq(from=xlow, to=xhigh, by=(xhigh-xlow)/5)
            graphics::axis(1, at=values, labels=round(values,2), las=1, 
                 font=plot_parameters[["font"]],
                 cex.axis=plot_parameters[["cex.axis"]], lwd.ticks=2) 
        }
        
    }
        
    #Mstat$concentration[Mstat$concentration==0] = 1
    
    graphics::points( Mstat$concentration, Mstat$value, pch = 19, 
            col=hsh_get( cond2color, as.character(Mstat$conditions_to_fit)), 
            cex=plot_parameters[["cex"]])
    if(show_legend){
        graphics::legend( x_legend, 0.45, x$unique_conditions, pch=19,
                col=hsh_get( cond2color, as.character(Mstat$conditions_to_fit)), 
                             bty="n", cex=0.75 )    
    }
    
    for(i in 1:dim(Mstat)[1]){
        px = Mstat$concentration[i]
        py = Mstat$value[i]
        pse = Mstat$sterr[i] * bar_multiple
        pcol = hsh_get( cond2color, as.character(Mstat$conditions_to_fit[i]) )
        bwdth = px/bandwidth_factor # varies if fit or not
        graphics::lines(c( px, px ), c( py-pse, py+pse ), lwd=1, col=pcol )
        graphics::lines(c( px-bwdth, px+bwdth ), c( py-pse, py-pse ), 
                        lwd=1, col=pcol )
        graphics::lines(c( px-bwdth, px+bwdth ), c( py+pse, py+pse ), 
                        lwd=1, col=pcol )        
    }
}

set_list_default = function(L, name, val){
    if( length( which( names(L)== name  ) )==0 )
        L[[name]] = val
    L
}
#' Plot interaction index with confidence intervals for observed effects 
#' at combination doses having observed effects.
#' 
#' Calls code published by Lee and Kong in Statistics in Biopharmaceutical 
#' Research 2012; please cite their work if you use this function.
#' 
#' @param ds dataset
#' @param sample_type sample type in ds
#' @param treatment_1 treatment in ds
#' @param treatment_2 treatment in ds
#' @param hour hour in ds. Default 0. 
#' @param log which scale should be logged; default is "y"
#' @param ... additional parameters to pass to plot function
#' @param alpha 1-alpha is the size of the confidence intervals, default 0.05
#' @return list ii: estimated interaction indices corresponding to the 
#' observations (c.d1, c.d2, E); ii.low, ii.up: estimated lower and upper CI
#' @references Lee & Kong Statistics in Biopharmaceutical Research 2012
#' @examples
#' dose_SCH=c(0.1, 0.5, 1, 2, 4)
#' eff_SCH = c(0.6701, 0.6289, 0.5577, 0.4550, 0.3755)
#' dose_4HPR = c(0.1, 0.5, 1, 2)
#' eff_4HPR = c(0.7666, 0.5833, 0.5706, 0.4934)
#' eff_comb = c(0.6539, 0.4919, 0.3551, 0.2341)
#' syn = data.frame( treatment_1 = rep("SCH66336", 13),
#'                   conc_1 = c( dose_SCH, rep(0, 4), dose_SCH[1:4]),
#'                   treatment_2 = rep("4-HPR", 13),
#'                   conc_2 = c( rep(0, 5), dose_4HPR, dose_4HPR ),
#'                   values = c(eff_SCH, eff_4HPR, eff_comb ), 
#'                   stringsAsFactors=FALSE )
#' ds_lk = create_synergy_dataset( sample_types = rep("sample_1", 13), 
#'                                 treatments_1 = syn$treatment_1,
#'                                 treatments_2 = syn$treatment_2,
#'                                 concentrations_1 = syn$conc_1,
#'                                 concentrations_2 = syn$conc_2,
#'                                 values = syn$values)
#' ii=plot_synergy_interaction_index(ds=ds_lk, sample_type = "sample_1", 
#'                               treatment_1 = "SCH66336", treatment_2="4-HPR", 
#'                               hour=0, 
#'                               xlab="Effect", ylab="Interaction Index" ,
#'                               main="Interaction Index")
#' @export
plot_synergy_interaction_index = function(ds, sample_type, 
                                          treatment_1, treatment_2, 
                                          hour, ..., log="y", alpha=0.05){
    
    # effects seen with combined treatment
    ds_type = is_dataset_valid(ds, treatments=c(treatment_1, treatment_2), 
                                hour=hour )
    pp = list(...)
    pp = set_list_default(pp, "xlim", c(0,1))
    pp = set_list_default(pp, "ylim", c(0.01, 10))
    pp = set_list_default(pp, "main", paste(sample_type, hour))
    pp = set_list_default(pp, "xlab", "Effect (Surviving fraction)")
    pp = set_list_default(pp, "ylab", "Interaction Index")
    pp = set_list_default(pp, "las", 1)
    pp = set_list_default(pp, "log", log)
    pp = set_list_default(pp, "pch", 19)
    
    CI.d=synergy_interaction_CI(ds, sample_type=sample_type, 
                                treatment_1=treatment_1, 
                                treatment_2=treatment_2,
                                hour=hour, alpha=alpha)
    y=CI.d$interaction_index
    x=CI.d$effects
    pp[["x"]] =  x[ CI.d$is_obs ]
    pp[["y"]] = y[ CI.d$is_obs ]
    do.call( graphics::plot, pp)
    graphics::lines( x, CI.d$cl_lower, lty=2, col="cornflowerblue" )
    graphics::lines( x, CI.d$cl_upper, lty=2, col="cornflowerblue" )
    graphics::lines( x, y, lty=1, col="black" )
    CI.d
}



#' Plot an image of the raw intensities for a plate
#' 
#' @param plate data frame in format matching that produced by 
#' \code{read_incucyte_from_excel}.
#' @param hour hour of the experiment to plot. Defaults to NA, plotting all 
#' values. If more than one time point has been stored for the data in plate, 
#' passing NA will result in an error. 
#' @param color_bounds values at which to draw coldest and hottest colors, 
#' defaults to c(0, 100)
#' @param color_palatte colors to use for cold-middle-hot, defaults to 
#' c("white", "#fdbb84","#e34a33")
#' @param main title to draw above plot, defaults to ""
#' @param cex point size, defaults to 1.5
#' @return none
#' @examples 
#' plate_1 = create_empty_plate( 6, hour=0, plate_id="plate_1")
#' plate_1[1,1:3] = c(100, 40, 10 )
#' plate_1[2,1:3] = c(90, 70, 30 )
#' plot_values_by_plate( plate_1, hour=0 )
#' # real-world example
#' pkg = "HTDoseResponseCurve"
#' fn_data = system.file("extdata", "sample_data_384.xlsx", package = pkg)
#' plate_data = read_plates_from_Incucyte_export( fn_data, "p1", 
#'                                                number_of_wells=384)
#' plot_values_by_plate(plate_data, hour=96)
#' @export
plot_values_by_plate = function( plate, hour=NA, color_bounds=c(0,100), 
                                 color_palatte=c("white", "#fdbb84","#e34a33"),
                                 main="", cex=1.5 ){
    if( is.na(hour) ){
        hour = unique(plate$hours)
        if( length(hour)>1 ){
            stop(paste("plate contains more than one time point; must specify",
                 "time point with the hour parameter") )
        }
    }
    if( sum( plate$hours==hour )==0 ){
        stop(paste("hour parameter", hour, "not present in plate"))
    }
    # plot raw value for a plate at the specified hour
    n_row = sum(plate$hours==hour)
    idx_max = dim(plate)[2] - 2
    n_col = idx_max
    graphics::plot( -1, -1, xlim=c(0, n_col+1), ylim=c(0, n_row+1), 
          axes=FALSE, xlab="", ylab="",yaxs="i", xaxs="i", main=main )
    abc = abc[1:n_row]
    abc_rev = abc[ length(abc):1]
    graphics::axis(1, 1:n_col, cex=0.5, las=2)
    graphics::axis(2, at=1:n_row, labels=abc_rev, las=1, cex=0.5)
    cmap = grDevices::colorRampPalette( colors=color_palatte )(color_bounds[2]-
                                                        color_bounds[1])
    vals = ceiling( data.matrix(plate[plate$hours==hour,1:idx_max] ) )
    colors = color_scale( vals, cmap, color_bounds = color_bounds )
    pches = matrix(19, ncol=n_col, nrow=n_row)
    pches[ is.na(colors) ] = 1
    colors[is.na(colors)] = "grey"
    xs = matrix( rep( 1:n_col, n_row ), ncol=n_col, nrow=n_row, byrow=TRUE)
    ys = matrix( rep( 1:n_row, n_col ), ncol=n_col, nrow=n_row)
    ys = n_row-ys+1
    graphics::points( as.numeric(xs), as.numeric(ys), col=colors, 
            pch=as.numeric(pches), cex=cex )
    graphics::box()
}

se = function(x){
    stats::sd(x, na.rm=TRUE)/sqrt(sum(!is.na(x))) 
}


#' Plot an time course of the raw intensities 
#' 
#' @param D experiment dataset with columns matching the output of 
#' \code{combine_data_and_maps}
#' @param sample_types which sample types to draw
#' @param treatments which treatments to draw
#' @param concentrations which concentrations to draw
#' @param plot_raw plot un-normalized raw values, defaults TRUE. If FALSE, 
#' plots value_normalized.
#' @param show_vehicle If TRUE and plot_raw is TRUE, plot un-normalized 
#' raw vehicle values. Defaults to FALSE. 
#' @param SEM_to_show If not zero, multiple of the SEM (standard error of 
#' the mean) to show at each point. Defaults to 0.
#' @param combine_raw_plates if TRUE, allow raw timecourse plots that include 
#' more than one plate. Normally combining more than one plate in a raw 
#' timecourse plot produces unexpected plot behavior because of plate-specific 
#' growth differences; that is not allowed unless combine_raw_plates is passed 
#' to specify that you know what you're doing. Defaults to FALSE.
#' @param ... standard parameters for \code{plot} function
#' @param cex.xaxis size coefficient for X axis, default 1
#' @param cex.yaxis size coefficient for Y axis, default 2
#' @param axis.font font value for axis, default 2 (bold)
#' @param summary_method summary method for points to plot in timecourse, one 
#' of (mean, median), defaults to mean
#' @return data frame reporting points plotted, sample types, and concentrations
#' @examples 
#' pkg = "HTDoseResponseCurve"
#' fn_map = system.file("extdata", "sample_data_384_platemap.txt",package=pkg)
#' fn_data = system.file("extdata", "sample_data_384.xlsx", package = pkg)
#' plate_map = read_platemap_from_Incucyte_XML( fn_map )
#' plate_data = read_plates_from_Incucyte_export( fn_data, "p1", 
#'                                                number_of_wells=384)
#' plate_data$hours = round(plate_data$hours)
#' ds = combine_data_and_map( plate_data, plate_map, negative_control = "DMSO" )
#' ds = ds[ds$treatment=="drug13" | ds$treatment=="DMSO",]
#' plot_timecourse( ds, sample_types=c("line_1", "line_2"), 
#'                      treatments="DMSO", concentrations=0)
#' @export
plot_timecourse = function( D, sample_types, treatments, 
                                concentrations, plot_raw=TRUE, 
                                show_vehicle=FALSE,
                                SEM_to_show = 0,
                                combine_raw_plates=FALSE,
                                ..., cex.xaxis=1, cex.yaxis=2, axis.font=2,
                                summary_method="mean"){
    if( summary_method != "mean" & summary_method != "median"){
        stop("summary_method parameter must be either mean or median")
    }
    for(i in 1:length(sample_types)){
        if( sum(D$sample_type==sample_types[i], na.rm=TRUE)==0 ){
            stop( paste("sample type", sample_types[i],
                        "in parameter sample_types not found in D") )
        }
    }
    for(i in 1:length(treatments)){
        if( sum(D$treatment==treatments[i], na.rm=TRUE)==0 ){
            stop( paste("treatment", treatments[i],
                        "in parameter treatments not found in D") )
        }
    }
    for(i in 1:length(concentrations)){
        if( sum(D$concentration==concentrations[i], na.rm=TRUE)==0 ){
            stop( paste("concentration",concentrations[i],
                        "in parameter concentrations not found in D") )
        }
    }
    ds_cur = D[ which(D$sample_type %in% sample_types &
                          D$treatment %in% treatments & 
                          D$concentration %in% concentrations),]
    
    if(dim(ds_cur)[1]==0){
        stop(paste("No matches for this combination of sample_types,",
                   "treatments, and concentrations"))
    }
    if( plot_raw & !combine_raw_plates & length(unique(D$plate_id))>1 ){
        stop(paste("if plot_raw is TRUE, you must restrict the dataset to a ",
                   "single plate unless combine_raw_plates is passed with",
                   "the value TRUE") )
    }
    if( show_vehicle ){
        # identify vehicle rows for each combination of sample_type+treatment
        veh = unique(ds_cur[,c("negative_control", "sample_type", "plate_id")])
        for( i in 1:dim(veh)[2] ){
            idx = which( D$treatment==veh$negative_control[i] & 
                         D$plate_id==veh$plate_id[i] & 
                         D$sample_type==veh$sample_type[i] )
            ds_cur = rbind(ds_cur, D[idx,] )   
        }
    }
    ds_cur = ds_cur[ !is.na( ds_cur$concentration ) & 
                         !is.na( ds_cur$treatment ) & 
                         !is.na( ds_cur$sample_type),]
    
    hours = sort( unique( ds_cur$hours ) )
    unique_concs = sort( unique(ds_cur$concentration) )
    unique_lines = sort( unique(sample_types) )
    unique_treatments = sort( unique (treatments ) )
    
    plot_parameters = list(...)
    if( length( which( names(plot_parameters)=="pch") ) > 0 ){
        if( length( plot_parameters[["pch"]]) != length(unique_lines) ){
            stop(paste("Passed pch parameter with",
                       length(plot_parameters[["pch"]]),
                       "elements but need to plot",length(unique_lines),
                       "unique lines"))
        }
        line2pch = hsh_from_vectors( as.character(unique_lines) , 
                                     plot_parameters[["pch"]] )
    }else{
        line2pch = hsh_from_vectors( as.character(unique_lines) , 
                                 INC_pches[1:length(unique_lines)] )
    }
    
    # eliminate NOTE generated by R CMD check 
    value=0; value_normalized=0; sample_type=0; concentration=0
    # If plotting raw data, plot one curve for each plate.
    Mstat = plyr::ddply( 
        .data=ds_cur, 
        .variables=c("sample_type", "hours", "concentration"), 
        summarize, 
        mu_raw = mean( value, na.rm=TRUE ),
        med_raw = stats::median( value, na.rm=TRUE ), 
        mu_norm = mean( value_normalized, na.rm=TRUE ),
        med_norm = stats::median( value_normalized, na.rm=TRUE ), 
        sem_raw = sd(value, na.rm=TRUE)/sqrt(sum(!is.na(value))) ,
        sem_norm= sd(value_normalized, na.rm=TRUE)/
                        sqrt(sum(!is.na(value_normalized))) ,
        sample_type=unique( sample_type ) ,
        concentration=unique( concentration )  )
    if( plot_raw ){
        if( summary_method=="mean"){
            values = Mstat$mu_raw
        }else{
            values = Mstat$med_raw
        }
        sems = Mstat$sem_raw
    }else{
        if( summary_method=="mean"){
            values = Mstat$mu_norm
        }else{
            values = Mstat$med_norm
        }
        sems = Mstat$sem_norm
    }
    
    conc2color = hsh_from_vectors( as.character(unique_concs) , 
                                   INC_colors_DRC[ 1:length(unique_concs) ] )
    Mstat = cbind( Mstat, 
                   color =hsh_get(conc2color,as.character(Mstat$concentration)), 
                   pch = hsh_get( line2pch, Mstat$sample_type),
                   stringsAsFactors=FALSE)
    
    plot_parameters[["xlim"]] = c(0, max(hours))
    plot_parameters[["axes"]] = FALSE
    plot_parameters[["pch"]] = Mstat$pch
    if( length( which( names(plot_parameters)=="col") ) == 0 ){
        plot_parameters[["col"]] = Mstat$color
    }
    plot_parameters[["x"]] = Mstat$hours
    plot_parameters[["y"]] = values
    if( length( which( names(plot_parameters)== "ylim"  ) )==0 ){
        plot_parameters[["ylim"]] = c(0, max(values, na.rm=TRUE) )
    }
    if( length( which( names(plot_parameters)== "ylab"  ) )==0 ){
        plot_parameters[["ylab"]] = ""
    }
    if( length( which( names(plot_parameters)== "xlab"  ) )==0 ){
        plot_parameters[["xlab"]] = ""
    }
    ylim = plot_parameters[["ylim"]]
    do.call( graphics::plot, plot_parameters )
    graphics::axis(2, round( c(ylim[1], (ylim[2])/2, ylim[2]), 1), las=1, 
                   cex.axis=cex.yaxis,font= axis.font )
    graphics::axis(1, round(hours,1), las=2, cex.axis=cex.xaxis, font=axis.font )
    graphics::box()
    if( SEM_to_show != 0 ){
        for( i in 1:length( plot_parameters[["x"]]  )){
            xx = plot_parameters[["x"]][i]
            yy = plot_parameters[["y"]][i] 
            lines( c(xx, xx), c( yy - (SEM_to_show*sems[i] ), 
                                 yy + (SEM_to_show*sems[i]) ), 
                   col=plot_parameters[["col"]][i] )
        }
    }
        
    
    unique(Mstat)
}


#' Plot a dose response curve summary values across one or more time points
#' 
#' Used to generate a line and dot plot of one or more AUC or EC50 values 
#' across a series of time points. All of the values in the fit_stats data 
#' frame will be plotted. Observations that do not meet the criteria specified 
#' in the alpha parameter are plotted using open values (e.g. pch=1). 
#' Observations where the minimum EC50 exceeds the obs_min parameter can be 
#' plotted using a user-specified value for the pch parameter; by default this 
#' is pch=13.
#' 
#' @param fit_stats data frame returned by call to \code{fit_statistics}
#' @param statistic statistic to plot, must be one of AUC, EC50
#' @param ... standard parameters for \code{plot} function
#' @param alpha cut-off to distinguish results with an ANOVA_P_value that is 
#' considered statistically significant, plotted as closed circles. Defaults to 
#' 1. 
#' @param obs_min cut-off to distinguish results where the minimum observed 
#' response across any dose is no greater than obs_min. Useful for restricting 
#' plots to those values that acheived a result such as a Surviving Fraction 
#' below 50 percent. Defaults to NA, which plots everything.
#' @return data frames with plotted points, colors. Useful for drawing a legend.
#' @examples 
#' pkg = "HTDoseResponseCurve"
#' fn_map = system.file("extdata","sample_data_384_platemap.txt",package=pkg)
#' fn_data = system.file("extdata", "sample_data_384.xlsx", package = pkg)
#' plate_map = read_platemap_from_Incucyte_XML( fn_map )
#' plate_data = read_plates_from_Incucyte_export( fn_data, "p1", 
#'                                                number_of_wells=384)
#' plate_data$hours = round(plate_data$hours)
#' ds = combine_data_and_map( plate_data, plate_map, negative_control = "DMSO" )
#' ds = ds[ds$treatment=="drug13",]
#' fits = fit_statistics(ds, fct = drc::LL.3() )
#' res=plot_fit_statistic( fits, "AUC", ylim=c(0, 6) ) 
#' @export
plot_fit_statistic = function( fit_stats, statistic, ..., alpha = 1, 
                               obs_min=NA){
    
    if( !( statistic == "AUC" | statistic=="EC50" ) ){
        stop("parameter statistic must be one of AUC, EC50" )
    }
    
    if( is.na(obs_min) ){
        obs_min = max(fit_stats$obs_min, na.rm=TRUE) + 1   
    }
    plot_par = list(...)
    if( length( which( names(plot_par)== "ylim"  ) )==0 ){
        plot_par[["ylim"]] = c(0, 10000)
    }
    if( length( which( names(plot_par)== "col"  ) )==0 ){
        colors = INC_colors_DRC    
    }else{
        colors = plot_par[["col"]]
    }
    if( length( which( names(plot_par)== "xlab"  ) )==0 ){
        plot_par[["xlab"]] = ""    
    }
    if( length( which( names(plot_par)== "cex"  ) )==0 ){
        plot_par[["cex"]] = 2
    }
    if( length( which( names(plot_par)== "ylab"  ) )==0 ){
        plot_par[["ylab"]] = ""    
    }
    plot_par[["x"]] = -100
    plot_par[["y"]] = -100
    if( length( which( names(plot_par)== "xlim"  ) )==0 ){
        plot_par[["xlim"]] = c(0, max(fit_stats$hour))
    }
    fit_stats = fit_stats[fit_stats$hour <= plot_par[["xlim"]][2],]
    plot_par[["axes"]] = FALSE
    ylim = plot_par[["ylim"]]
    
    if( length( which( names(plot_par)== "xlim"  ) )==0 ){
        plot_par[["xlim"]] = c(0, max(fit_stats$hour, na.rm=TRUE) )
    }
    
    do.call( graphics::plot, plot_par )
    ylims = 1:ylim[2]
    line_y = seq(from=ylim[1], to=ylim[2], by=(ylim[2]-ylim[1])/5 )
    for(i in 1:length(line_y)){
        graphics::abline( line_y[i], 0, col="lightgray")   
    }
    
    if( length( which( names(plot_par)== "las"  ) )==0 ){
        plot_par[["las"]] = 1
    }
    
    graphics::axis(2, line_y, las=plot_par[["las"]])
    graphics::axis(1, sort(unique(fit_stats$hour)), las=plot_par[["las"]])
    ctr=1
    sample_types = sort(unique(fit_stats$sample_type))
    treatments = sort(unique(fit_stats$treatment))
    for( i in 1:length( sample_types ) ){
        cur_samp = sample_types[i]
        for( j in 1:length(treatments) ){
            cur_treat = treatments[j]
            fits_cur = fit_stats$sample_type==cur_samp & 
                fit_stats$treatment==cur_treat
            if( statistic=="AUC" ){
                y=fit_stats$AUC[ fits_cur ]
            }else if( statistic=="EC50" ){
                y=fit_stats$coef_EC50[ fits_cur ]
            }else{
                stop(paste("cannot plot statistic",statistic))
            }
            xs = unique(fit_stats$hour[ fits_cur ] )
            pvals = fit_stats$ANOVA_P_value[ fits_cur ]
            min_observations = fit_stats$obs_min[ fits_cur ]
            if( length( y ) >1 ){
                for(k in 2:length(y)){
                    if( !is.na( pvals[k-1] ) & !is.na( pvals[k] ) & 
                                pvals[k-1] <= alpha & pvals[k] <= alpha ){
                        graphics::lines( c(xs[k], xs[k-1] ), 
                               c( y[k],  y[k-1] ), col=colors[ctr] )
                    }
                }
            }
            
            pches = rep(19, length(pvals))
            pches[ pvals > alpha ] = 13    
            pches[ pvals <= alpha & min_observations > obs_min ] = 18 
            
            graphics::points( xs, y, col=colors[ctr] , pch=pches, 
                    cex=plot_par[["cex"]])
            for(i in 1:length(y)){
                if( is.infinite( y[i] ) ){
                    graphics::text( xs[i], ylim[2], "Inf", cex=1, font=2, 
                                    col=colors[ctr])
                }
            }
            cur_res = data.frame( sample_type=rep(cur_samp, length(y)),
                                  treatment=rep(cur_treat, length(y)),
                                  AUC=y, hour=xs, 
                                  color=rep(colors[ctr], length(y) ), 
                                  stringsAsFactors=FALSE)
            if( ctr==1 ){
                res = cur_res    
            }
            else{
                res = rbind(res, cur_res)   
            }
            ctr = ctr+1
        }
    }
    res
}



#' Standard boxplot with labeled outliers
#' 
#' Any point outside of the whiskers will be labeled. Passes all standard 
#' \code{boxplot()} parameters through to the R \code{boxplot()} function.
#' 
#' @param M one-dimensional matrix of values with row names for labels
#' @param ... standard commands for \code{boxplot()} function
#' @return standard output of \code{boxplot()} function
#' @examples
#' M = matrix( c(3,4,7,5,6,10), nrow=6, ncol=1 )
#' dimnames(M)[[1]] = paste("drug",1:6)
#' boxplot_label_outliers(M[,1])
#' @export
boxplot_label_outliers = function( M, ... ){
    plot_parameters = list(...)
    if( length( which( names(plot_parameters)== "ylim"  ) )==0 )
        plot_parameters[["ylim"]] = c(0, ceiling(max(M, na.rm=TRUE) ) )
    if( length( which( names(plot_parameters)== "xlim" ) )== 0 )
        plot_parameters[["xlim"]] = c(0.25, 1.25)
    y_max = plot_parameters[["ylim"]][2]
    y_min = plot_parameters[["ylim"]][1]
    if( length( which( names(plot_parameters)== "cex"))==0 )
        plot_parameters[["cex"]] = 1
    
    plot_parameters[["x"]] = M 
    b=do.call( graphics::boxplot, plot_parameters )
    if( length(b$out)>0 ){
        outs = b$out
        outs = sort(outs, decreasing=TRUE)
        xs = rep(1, length(M))
        not_out = which( !( round(M,3) %in% round(outs,3)) )
        xs[not_out] = jitter( xs[not_out], 2 )
        graphics::points( xs, M, pch=19, cex=0.5 )
        interval = (y_max-y_min) / (length(b$out)+1)
        y_cur = y_max-interval
        for( i in 1:length(b$out)){
            graphics::text( 0.7, y_cur, names(outs)[i], adj=1, 
                  cex=plot_parameters[["cex"]] )
            graphics::lines( c(0.7,1), c(y_cur, outs[i]), col="#00000033" ) 
            y_cur = y_cur-interval
        }
    }
    b
}



#' plot FA vs. CI Chou synergy median effect plot
#' @param CS chou statistics
#' @param ... standard parameters for \code{plot} function
#' @return list of linear models fit
#' @examples 
#' dose_SCH=c(0.1, 0.5, 1, 2, 4)
#' eff_SCH = c(0.6701, 0.6289, 0.5577, 0.4550, 0.3755)
#' dose_4HPR = c(0.1, 0.5, 1, 2)
#' eff_4HPR = c(0.7666, 0.5833, 0.5706, 0.4934)
#' eff_comb = c(0.6539, 0.4919, 0.3551, 0.2341)
#' syn = data.frame( 
#'     treatment_1 = rep("SCH66336", 13),
#'     conc_1 = c( dose_SCH, rep(0, 4), dose_SCH[1:4]),
#'     treatment_2 = rep("4-HPR", 13),
#'     conc_2 = c( rep(0, 5), dose_4HPR, dose_4HPR ),
#'     values = c(eff_SCH, eff_4HPR, eff_comb ),
#'     stringsAsFactors=FALSE )
#' ds_lk = create_synergy_dataset( sample_types = rep("sample_1", 13), 
#'                                 treatments_1 = syn$treatment_1,
#'                                 treatments_2 = syn$treatment_2,
#'                                 concentrations_1 = syn$conc_1,
#'                                 concentrations_2 = syn$conc_2,
#'                                 values = syn$values)
#' CS_lk = chou_synergy( ds_lk, sample_type = "sample_1", hour=0,
#'                       treatment_1 = "SCH66336", treatment_2="4-HPR",
#'                       fct=drc::LL.2(), summary_method="mean" )
#' me=plot_synergy_median_effect( CS_lk, main="Median Effect" )
#' @export
plot_synergy_median_effect=function( CS, ... ){   
    cs1 = CS$treatment_1
    cs2 = CS$treatment_2
    csc = CS$treatment_12
    cs1$x_m[cs1$D==0] = NA
    cs1$y_m[cs1$D==0] = NA
    cs2$x_m[cs2$D==0] = NA
    cs2$y_m[cs2$D==0] = NA
    csc$x_m[csc$D==0] = NA
    csc$y_m[csc$D==0] = NA
    ymax = ceiling( max( c( abs(cs1$y_m),abs(cs2$y_m),abs(csc$y_m)),na.rm=TRUE))
    xmax = ceiling( max( c( abs(cs1$x_m),abs(cs2$x_m),abs(csc$x_m)),na.rm=TRUE))

    plot_parameters = list(...)
    if( length( which( names(plot_parameters) == "ylim"))==0 )
        plot_parameters[["ylim"]] = c(-1*ymax, ymax)
    if( length( which( names(plot_parameters) == "xlim"))==0 )
        plot_parameters[["xlim"]] = c(-1*xmax, xmax)
    if( length( which( names(plot_parameters) == "ylab"))==0 )
        plot_parameters[["ylab"]] = "Log(Fa/Fu)"
    if( length( which( names(plot_parameters) == "xlab"))==0 )
        plot_parameters[["xlab"]] = "Log(dose)"
    if( length( which( names(plot_parameters) == "las"))==0 )
        plot_parameters[["las"]] = 1
    if( length( which( names(plot_parameters) == "col"))==0 )
        plot_parameters[["col"]] = c("black", "red", "cornflowerblue")
    if( length( which( names(plot_parameters) == "lwd"))==0 )
        plot_parameters[["lwd"]] = 1
    if( length( which( names(plot_parameters) == "pch"))==0 )
        plot_parameters[["pch"]] = 19
    if(length(plot_parameters[["col"]]) != 3 ){
        stop("must pass three colors in col argument")
    }
    colors = plot_parameters[["col"]]
    pches = plot_parameters[["pch"]]
    if( length(pches)==1 )
        pches = rep( pches[1], 3)
    plot_parameters[["x"]] = cs1$x_m
    plot_parameters[["y"]] = cs1$y_m
    plot_parameters[["col"]] = colors[1]
    plot_parameters[["pch"]] = pches[1]
    do.call( graphics::plot, plot_parameters)
    graphics::lines(c(0,0), c(-100,100))
    graphics::abline(0,0)
    graphics::points( cs2$x_m, cs2$y_m, pch=pches[2], col=colors[2])
    graphics::points( csc$x_m, csc$y_m, pch=pches[3], col=colors[3])
    LM_1 = stats::lm(cs1$y_m~cs1$x_m)
    LM_2 = stats::lm(cs2$y_m~cs2$x_m)
    LM_12 = stats::lm(csc$y_m~csc$x_m)
    graphics::abline( LM_1, col=colors[1], lwd=plot_parameters[["lwd"]] )
    graphics::abline( LM_2, col=colors[2], lwd=plot_parameters[["lwd"]]  )    
    graphics::abline( LM_12, col=colors[3], lwd=plot_parameters[["lwd"]]  )
    list( LM_1=LM_1, LM_2=LM_2, LM_12=LM_12)
}

#' plot FA vs. CI Chou synergy plot
#' @param CS chou statistics
#' @param ... standard parameters for \code{plot} function
#' @param show_horizontal show horizontal stripes
#' @examples 
#' dose_SCH=c(0.1, 0.5, 1, 2, 4)
#' eff_SCH = c(0.6701, 0.6289, 0.5577, 0.4550, 0.3755)
#' dose_4HPR = c(0.1, 0.5, 1, 2)
#' eff_4HPR = c(0.7666, 0.5833, 0.5706, 0.4934)
#' eff_comb = c(0.6539, 0.4919, 0.3551, 0.2341)
#' syn = data.frame( treatment_1 = rep("SCH66336", 13),
#'                   conc_1 = c( dose_SCH, rep(0, 4), dose_SCH[1:4]),
#'                   treatment_2 = rep("4-HPR", 13),
#'                   conc_2 = c( rep(0, 5), dose_4HPR, dose_4HPR ),
#'                   values = c(eff_SCH, eff_4HPR, eff_comb ),
#'                   stringsAsFactors=FALSE )
#' ds_lk = create_synergy_dataset( sample_types = rep("sample_1", 13), 
#'                                 treatments_1 = syn$treatment_1,
#'                                 treatments_2 = syn$treatment_2,
#'                                 concentrations_1 = syn$conc_1,
#'                                 concentrations_2 = syn$conc_2,
#'                                 values = syn$values)
#' CS_lk = chou_synergy( ds_lk, sample_type = "sample_1", hour=0,
#'                       treatment_1 = "SCH66336", treatment_2="4-HPR",
#'                       fct=drc::LL.2(), summary_method="mean" )
#' plot_chou_synergy_Fa_CI(CS_lk)
#' @return none
#' @export
plot_chou_synergy_Fa_CI = function( CS, ..., show_horizontal=TRUE  ){
    plot_parameters = list(...)
    if( length( which( names(plot_parameters) == "ylab"))==0 )
        plot_parameters[["ylab"]] = "CI"
    if( length( which( names(plot_parameters) == "xlab"))==0 )
        plot_parameters[["xlab"]] = "Fa"
    if( length( which( names(plot_parameters) == "col"))==0 )
        plot_parameters[["col"]] = "red"
    if( length( which( names(plot_parameters) == "ylab"))==0 )
        plot_parameters[["ylab"]] = "CI"
    if( length( which( names(plot_parameters) == "ylim"))==0 )
        plot_parameters[["ylim"]] = c(0,2)
    if( length( which( names(plot_parameters) == "xlim"))==0 )
        plot_parameters[["xlim"]] = c(0,1)
    if( length( which( names(plot_parameters) == "pch"))==0 )
        plot_parameters[["pch"]] = 19
    plot_parameters[["x"]] = CS$CI$Fa
    plot_parameters[["y"]] = CS$CI$CI
    do.call( graphics::plot, plot_parameters )
    if(show_horizontal){
        for(i in seq(0, 2, 0.1)){
            graphics::abline(i,0, col="lightgray")
        }
        graphics::abline(1,0, lwd=2)
    }
    graphics::points( CS$CI$Fa, CS$CI$CI, pch=plot_parameters[["pch"]], 
            col=plot_parameters[["col"]] )
}

#' Plot an time course of the raw intensities 
#' 
#' @param ds synergy experiment dataset with columns matching the output of 
#' \code{combine_data_and_maps}
#' @param sample_type which sample type to draw; only one can be specified
#' @param treatment_a which of two synergy treatments is labeled treatment A
#' @param treatment_b which of two synergy treatments is labeled treatment B 
#' @param hour which hour to use from dataset ds, defaults to 0
#' @param normalize plot isobologram normalized so that the EC50 for each axis
#' is normalized to value 1, with other concentrations relative to this value
#' @param xmax highest value for X axis (treatment A)
#' @param ymax highest value for Y axis (treatment B)
#' @return list( fit_a, fit_b, IC50_a, IC50_b, iso) where fit_a and fit_b are 
#' the fit of treatments A and B alone, IC50_a and IC50_b, and isobologram, a 
#' data frame of the points X and Y that were plotted
#' @export 
plot_isobologram = function( ds, sample_type, treatment_a, treatment_b, hour=0,
                             normalize=TRUE, xmax=NA, ymax=NA){
    if( length(sample_type) != 1 ){
        stop("must call with a single sample type")   
    }
    TA = treatment_a
    TB = treatment_b
    ds_a = standardize_treatment_assigments( ds, TA, TB )
    ds_b = standardize_treatment_assigments( ds, TB, TA )
    ds_a_only = ds_a[ds_a$treatment != TB & ds_a$treatment_2 != TB,]
    ds_b_only = ds_b[ds_b$treatment != TA & ds_b$treatment_2 != TA,]
    fit_a = fit_DRC(ds_a_only, sample_types = c(sample_type), treatments = TA, 
                    fct=LL.3(), hour = hour)
    fit_b = fit_DRC(ds_b_only, sample_types = c(sample_type), treatments = TB, 
                    fct=LL.3(), hour = hour)
    IC50_a = fit_a$fit_stats$coef_EC50
    IC50_b = fit_b$fit_stats$coef_EC50
    if( normalize ){
        xmax=2
        ymax=2
        plot( 1, 0, xlim=c(0, xmax), ylim=c(0, ymax), xaxs="i", yaxs="i", 
              pch=19, xlab=treatment_a, ylab=treatment_b)
        graphics::points( 0, 1, pch=19)
        graphics::lines( c(0,1), c(1,0), lwd=2)
    }else{
        if( is.na(xmax) ){
            xmax = max( ds_a$concentration[ds$treatment==TA], na.rm=TRUE )
        }
        if( is.na(ymax) ){
            ymax = max( ds_b$concentration[ds$treatment==TB], na.rm=TRUE )   
        }
        plot( IC50_a, 0, xlim=c(0, xmax), ylim=c(0, ymax), xaxs="i", yaxs="i", 
              pch=19, xlab=treatment_a, ylab=treatment_b)
        graphics::points( 0, IC50_b, pch=19)
        graphics::lines( c(0, IC50_a), c(IC50_b, 0), lwd=2 )
    }
    # fit concentration for each concentration_2, 
    # estimate IC50, plot conc_IC50, conc_2
    concs_a = sort(unique(ds_a$concentration[!ds_a$is_negative_control]))
    concs_b = sort(unique(ds_b$concentration[!ds_b$is_negative_control]))
    iso_x = c()
    iso_y = c()
    for(i in 1:length(concs_b)){
        ds_iso=convert_dataset_for_synergy_DRC( ds_a[ds_a$hours==hour,],
                                                sample_type=sample_type,
                                                treatment_a = TA, treatment_b = TB,
                                                concentrations_a = concs_a,
                                                concentrations_b= concs_b[i])
        fit_par = fit_DRC(ds_iso, 
                          sample_types = unique(ds_iso$sample_type[!ds_iso$is_negative_control]),
                          treatments = TA,
                          fct=LL.4(), hour = hour)
        iso_x = c( iso_x, fit_par$fit_stats$coef_EC50[2] )
        iso_y = c( iso_y, concs_b[i] )
    }
    for(i in 1:length(concs_a)){
        ds_iso=convert_dataset_for_synergy_DRC( ds_b[ds_b$hours==hour,],
                                                sample_type=sample_type,
                                                treatment_a = TB, treatment_b = TA,
                                                concentrations_a = concs_b,
                                                concentrations_b= concs_a[i])
        fit_par = fit_DRC(ds_iso, 
                          sample_types = unique(ds_iso$sample_type[!ds_iso$is_negative_control]),
                          treatments = TB,
                          fct=LL.4(), hour = hour)
        iso_x = c( iso_x, fit_par$fit_stats$coef_EC50[2] )
        iso_y = c( iso_y, concs_a[i] )
    }
    
    if( normalize ){
        iso_x = iso_x / IC50_a
        iso_y = iso_y / IC50_b
    }
    
    points( iso_x, iso_y )
    iso = data.frame( x=iso_x, y=iso_y )
    list( fit_a, fit_b, IC50_a, IC50_b, isobologram=iso)
}

#' Plot additive vs. synergy effect for a given concentration of one drug
#' 
#' @param ds synergy experiment dataset with columns matching the output of 
#' \code{combine_data_and_maps}
#' @param treatment_for_DRC the treatment to vary across different doses
#' @param concs_for_DRC concentrations of treatment_for_DRC to show on X axis
#' @param treatment_subgroup the treatment to hold constant
#' @param conc_subgroup concentration of subgroup treatment
#' @param hour which hour to use from dataset ds, defaults to 0
#' @param sample_type sample type to draw; only one can be specified
#' @param ... accept standard arguments to plot()
#' @export 
plot_additive_vs_synergy_effect = function( ds, 
                                            treatment_for_DRC, concs_for_DRC, 
                                            treatment_subgroup, conc_subgroup, 
                                            hour, sample_type, ... ){
    # standardize treatment assignments and restrict to relevant time/sample
    ds_std = ds[ds$hours==hour & ds$sample_type==sample_type,]
    ds_std = standardize_treatment_assigments(ds_std, treatment_for_DRC, 
                                                          treatment_subgroup )
    ds_std=ds_std[ c(ds_std$concentration %in% concs_for_DRC & 
                     ds_std$concentration_2 %in% conc_subgroup )|
                     ds_std$is_negative_control | ds_std$is_negative_control_2,]
    
    # summarize replicates
    M=plyr::ddply( ds_std, c("treatment", "concentration", 
                             "treatment_2", "concentration_2",
                             "is_negative_control", "is_negative_control_2"), 
                   function(x){ data.frame( 
                              mu=round(mean(x$value_normalized, na.rm=TRUE),2) ) 
                   } )
    
    # single treatments, putting vehicle first
    Ma = M[ M$is_negative_control_2,]
    Ma = Ma[ order( !Ma$is_negative_control, Ma$concentration ),]
    Mb = M[ M$is_negative_control,]
    Mb = Mb[ order( !Mb$is_negative_control_2, Mb$concentration_2 ),]
    Mcomb = M[ (M$is_negative_control & M$is_negative_control_2) | 
               (!M$is_negative_control & !M$is_negative_control_2), ]
    Mcomb = Mcomb[ order( !Mcomb$is_negative_control_2, Mcomb$concentration ),]
    effect_both = 1 - Mcomb$mu
    effect_DRC = 1 - Ma$mu
    effect_subgroup = 1 - Mb$mu
    
    idx_subgroup = which( Mb$concentration_2 == conc_subgroup )
    
    plot_parameters = list(...)

    if( length( which( names(plot_parameters) == "main"))==0 )
        plot_parameters[["main"]] = paste(treatment_subgroup, conc_subgroup)    
    if( length( which( names(plot_parameters) == "ylab"))==0 )
        plot_parameters[["ylab"]] = 'effect'
    if( length( which( names(plot_parameters) == "xlab"))==0 )
        plot_parameters[["xlab"]] = treatment_for_DRC
    if( length( which( names(plot_parameters) == "col"))==0 ){
        plot_parameters[["col"]] = c("black","black","black")
    }
    if( length( which( names(plot_parameters) == "cex"))==0 ){
        plot_parameters[["cex"]] = 1
    }
    else if( length( plot_parameters[["col"]] )==1 ){
        color=plot_parameters[["col"]]
        plot_parameters[["col"]] = c(color, color, color)
    }else if( length( plot_parameters[["col"]] ) == 3 ){
            plot_parameters[["col"]] = plot_parameters[["col"]]
    }else{
        stop("Must pass either zero, one, or three colors")
    }
    if( length( which( names(plot_parameters) == "ylim"))==0 )
        plot_parameters[["ylim"]] = c(0,1)
    plot_parameters[["yaxs"]] = "i"
    plot_parameters[["x"]]=effect_DRC
    plot_parameters[["axes"]]=FALSE
    b=do.call( graphics::plot, plot_parameters )
    graphics::points( rep( effect_subgroup[idx_subgroup], length(effect_DRC)), 
                      pch=2, col=plot_parameters[["col"]][2], 
                      cex=plot_parameters[["cex"]] )
    graphics::points( effect_both, pch=19,col=plot_parameters[["col"]][3],
                      cex=plot_parameters[["cex"]] )
    axis(2, seq(from=0, to=1, by=.20), las=1)
    axis(1, at=1:(1+length(concs_for_DRC)), labels=c(0, concs_for_DRC), las=2 )
}


#' Create complete report for a single treatment
#' 
#' @param dataset dataset with columns matching the output of 
#' \code{combine_data_and_maps}
#' @param fit_summary output of a call to \code{fit_statistics}
#' @param treatment_for_DRC the treatment to display
#' @param sample_type sample types to show; if NA, show all in arbitrary order
#' @param times_DRC time points at which to plot dose response curves (maximum
#' of four, must be present in the dataset. If NA, no curves are plotted.
#' @param ylim_raw lower and upper y-axis limit for raw DMSO curves
#' @export 
#' 
plot_treatment_summary = function( dataset, fit_summary, 
                                   treatment, sample_types=NA,
                                   times_DRC=NA,
                                   ylim_raw=NA ){
    MAX_SF50 = 20e6
    if( is.na(ylim_raw)[1] ){
        ylim_raw = c(0, ceiling(max(dataset$value, na.rm=TRUE)) )
    }
    if( length(ylim_raw) != 2 ){
        stop("ylim_raw must be a vector of length 2")   
    }
    if( is.na(sample_types[1]) ){
        sample_types = sort( unique( dataset$sample_type) )
    }
    if( sum(fit_summary$treatment==treatment )==0 ){
        stop(paste("treatment",treatment,"not in fit_summary"))        
    }

    for(i in 1:length(sample_types)){
        if( sum(fit_summary$sample_type==sample_types[i])==0 ){
            stop(paste("sample_type",sample_types[i],"not in fit_summary"))        
        }
    }
    n_treat = is_dataset_valid(dataset, treatment=treatment, 
                               sample_types=sample_types)
    if( !is.na(times_DRC)[1]){
        if( length( times_DRC )>4 ){
            stop("Can only highlight four times for DRC plot.")   
        }
        for(i in 1:length(times_DRC)){
            if( sum(fit_summary$hour==times_DRC[i])==0 ){
                stop(paste("time",times_DRC[i],"not in fit_summary"))        
            }
        }
    }
    vehicle = unique( dataset$negative_control[dataset$treatment==treatment ] )
    if( sum(dataset$treatment==vehicle)==0){
        stop( "dataset does not contain vehicle measurements" )
    }
    ds = dataset[dataset$sample_type %in% sample_types &
                 dataset$treatment %in% c(treatment, vehicle),]
    plates = sort(unique(ds$plate_id[ds$treatment==treatment]))
    if( length(plates)>4 ){
        warning("Cannot currently plot raw data for more than 4 plates")   
    }
    if( length(plates)==1 ){ n_blanks = 3 }
    if( length(plates)==2 ){ n_blanks = 2 }
    if( length(plates)==3 ){ n_blanks = 1 }
    
    # Axis should have all hours in the experiment
    # individual plates may be missing timepoints, so need to re-calculate the
    # hours used to plot values for each metric
    AXIS_hours = sort(unique(round(ds$hours[ds$treatment==treatment],1)))
    if( is.na( times_DRC ) ){
        layout(matrix( 1:8, byrow=TRUE, ncol=4, nrow=2 ))
    }else{
        layout(matrix( 1:12, byrow=TRUE, ncol=4, nrow=3 ))
    }
    par(mar=c(3,5,2,0))
    for( i in 1:length(plates) ){
        conc_this_plate = min(unique(ds$concentration[ds$plate_id==plates[i] & 
                                                      ds$treatment==vehicle]))
        samp_this_plate = sort(unique(ds$sample_type[ds$plate_id==plates[i] & 
                                                      ds$treatment==vehicle]))
        #ymax = ceiling( ds$value[ds$plate_id==plates[i] & ds$treatment==treatment )
        tc=plot_timecourse(ds[ds$plate_id==plates[i],],
                          sample_types = samp_this_plate ,
                          treatments = vehicle,
                          plot_raw = TRUE,
                          show_vehicle = TRUE,
                          concentrations = conc_this_plate,
                          cex.yaxis = 1, 
                          ylim=ylim_raw,
                          main=paste("raw", vehicle, plates[i] ))
        labels = unique(tc[,c("sample_type", "pch")])$sample_type
        pches = unique(tc[,c("sample_type", "pch")])$pch
        legend(0, ylim_raw[2], labels, cex=0.75, pch=pches, bty = "n")
    }
    
    for(i in 1:n_blanks){
        plot( -1000, -1000, axes=FALSE, xlab="", ylab="", xlim=c(0,1),
              ylim=c(0, 1), main="")
    }
    
    
    ymax = 1.4
    if( max(fit_summary$AUC[fit_summary$treatment==treatment], na.rm=TRUE)>1.5){ 
        ymax = 2
    }
    plot( -1000, -1000, axes=FALSE, xlab="", ylab="AUC", 
          xlim=c(0, max(AXIS_hours)), ylim=c(0, ymax), 
          main=paste( "AUC", treatment) )
    
    graphics::axis(2, 0:ceiling(ymax), las=2, cex.axis=1, font=2 )
    graphics::axis(1, AXIS_hours, las=2, cex.axis=1, font=2 )
    for(i in 1:length(sample_types)){
        idx = fit_summary$sample_type==sample_types[i] & 
              fit_summary$treatment==treatment
        vals = fit_summary$AUC[ idx ]
        hours = round( fit_summary$hour[ idx ], 1)
        points( hours, vals, col=INC_colors_DRC[i], pch=19 )
    }
    legend(0, ymax, sample_types, cex=0.75, pch=19, 
           col=INC_colors_DRC[1:length(sample_types)], bty = "n")
    box()

    ymax = 1.4
    MAX_EC50 = 20e6
    
    EC50 = fit_summary$coef_EC50
    EC50[ is.infinite(EC50)] = NA
    EC50[ EC50 > MAX_EC50 ] = MAX_EC50
    ymax = ceiling( max( log10( EC50[fit_summary$treatment==treatment]), na.rm=TRUE) )
    EC50 = log10(EC50)
    plot( -1000, -1000, axes=FALSE, xlab="", ylab="log EC50", 
          xlim=c(0, max(AXIS_hours)), ylim=c(0, ymax), 
          main=paste( "EC50", treatment))
    labels = c()
    graphics::axis(2, las=2, cex.axis=1, font=2,
                at = 0:ceiling(ymax),
                labels=parse( text=paste("10^",0:ceiling(ymax),sep="")) )
    graphics::axis(1, AXIS_hours, las=2, cex.axis=1, font=2 )
    for(i in 1:length(sample_types)){
        idx = fit_summary$sample_type==sample_types[i] & 
              fit_summary$treatment==treatment
        vals = EC50[ idx ]
        hours = round( fit_summary$hour[ idx ], 1)
        points( hours, vals, col=INC_colors_DRC[i], pch=19 )
        if(length(!is.na(vals))>0){
            points( round(hours,1), vals, col=INC_colors_DRC[i], pch=19 )
        }
    }
    legend(0, ymax, sample_types, cex=0.75, pch=19, 
           col=INC_colors_DRC[1:length(sample_types)], bty = "n")
    box()
    
    ymax = 1.4
    SF50 = fit_summary$SF50
    SF50[ is.infinite(SF50)] = NA
    SF50[ SF50 > MAX_SF50 ] = MAX_SF50
    ymax = ceiling( max( log10( SF50[fit_summary$treatment==treatment]), na.rm=TRUE) )
    if( is.na(ymax) | is.infinite(ymax) )
        ymax=6
    SF50 = log10(SF50)
    plot( -1000, -1000, axes=FALSE, xlab="", ylab="log SF50", 
          xlim=c(0, max(hours)), ylim=c(0, ymax), 
          main=paste( "SF50", treatment))
    graphics::axis(2, las=2, cex.axis=1, font=2,
                   at = 0:ceiling(ymax),
                   labels=parse( text=paste("10^",0:ceiling(ymax),sep="")) )
    graphics::axis(1, AXIS_hours, las=2, cex.axis=1, font=2 )
    for(i in 1:length(sample_types)){
        idx = fit_summary$sample_type==sample_types[i] & 
            fit_summary$treatment==treatment
        vals = SF50[ idx ]
        hours = round( fit_summary$hour[ idx ], 1)
        points( hours, vals, col=INC_colors_DRC[i], pch=19 )
        if( length(!is.na(vals))>0){
            points( round(hours,1), vals, col=INC_colors_DRC[i], pch=19 )
        }
    }
    legend(0, ymax, sample_types, cex=0.75, pch=19, 
           col=INC_colors_DRC[1:length(sample_types)], bty = "n")
    box()
    
    # All ANOVA p values will be same at a given {time point, treatment}
    idx = fit_summary$sample_type==sample_types[1] & 
          fit_summary$treatment==treatment
    PVAL = -1 * log10( fit_summary$ANOVA_P_value[idx] )
    PVAL[ is.infinite(PVAL)] = NA
    vals = PVAL
    ymax = ceiling( max(vals, na.rm=TRUE) )
    hours = round( fit_summary$hour[ idx ], 1)
    plot( hours,
          vals,
          axes = FALSE, xlab="", ylab="-log(P)", 
          xlim = c( 0, max(AXIS_hours) ), ylim=c( 0, ymax ), 
          main = paste( "ANOVA", treatment), pch=19 )
    graphics::axis(2, las=2, cex.axis=1, font=2,
                   at = 0:ceiling(ymax),
                   labels=parse( text=paste("10^",0:ceiling(ymax),sep="")) )
    graphics::axis(1, AXIS_hours, las=2, cex.axis=1, font=2 )
    box()
    
    if( !is.na(times_DRC)){
        for(h in 1:length(times_DRC)){
            fit=fit_DRC( ds, 
                     sample_types=sample_types,
                     treatments = treatment, 
                     hour=times_DRC[h], 
                     fct=drc::LL.3() )
            fpc=HT_fit_plot_colors( fit ) 
            plot( fit, 
              xlim=c(0, 1e6), ylim=c(0, 1.2),
              lwd=3, 
              main=paste(treatment, times_DRC[h], "hr"), 
              ylab="surviving fraction", 
              xlab="nM")
            legend( 1, 0.35, 
                legend = get.split.col(fpc$condition, "_|_", first = TRUE),
                pch=19, col=fpc$col, bty="n", cex=0.75)
            box()
        }
        if( length(times_DRC)<4 ){
            for(i in 1 : (4-length(times_DRC)) ){
                plot( -1000, -1000, axes=FALSE, xlab="", ylab="", xlim=c(0,1),
                      ylim=c(0, 1), main="")
            }
        }
    }
    
}