<jupyter_start><jupyter_text>Preamble<jupyter_code>libraries = c("dplyr","magrittr","tidyr","ggplot2","gridExtra","RColorBrewer") #,"zoo","directlabels")
for(x in libraries) {
    library(x,character.only=TRUE,warn.conflicts=FALSE) }

# options(jupyter.plot_mimetypes = "image/svg+xml") 
clrs = brewer.pal(8,"Set1")

windowsFonts(Times = windowsFont("Times New Roman"))
base_sz = 12 # base_size parameter

'%&%' = function(x,y)paste0(x,y)
'%!in%' = function(x,y)!('%in%'(x,y))
    
# Initialization of array for recorded plots
plot_point_sizes = c(); nm = c(); plot_point_sizes = list()<jupyter_output><empty_output><jupyter_text>Sensitivity analysis for $\mu$ and $\bar\mu$<jupyter_code>"../results/final-sensitivity_mu.csv" %>% 
    read.csv %>%
    select(-matches("err|sw_end|univ_point")) %>%
    arrange(T,Tbar) -> df
#     group_by(T) %>% mutate(y = 1:n()) %>% ungroup %>%
#     group_by(Tbar) %>% mutate(x = 1:n()) %>% ungroup 
df %>% head<jupyter_output><empty_output><jupyter_text>Main figure for $c=4\%$<jupyter_code>cs = c(4.2,3.5)
options(repr.plot.width=cs[1],repr.plot.height=cs[2])

outcome_threshold = 10
nbin=20; w=3; h=2.2
cl_rc = brewer.pal(8,"Greys")[2]
cl_pt = "green" #brewer.pal(8,"Greens")[3]
vjst = .4
hjst = .5
df %>% filter(index==1) %>%
    mutate(outcome_cut=if_else(outcome>=outcome_threshold,outcome_threshold,outcome)) %>%
    ggplot(aes(x=T,y=Tbar,fill=outcome_cut)) + 
        geom_raster(interpolate=T) +
        geom_abline(slope=1,size=.2,linetype="dashed") + 
        stat_contour(aes(z=outcome),breaks=c(.2,.5,2,5), color=brewer.pal(8,"Greys")[1], size=.9) +
        stat_contour(aes(z=outcome),breaks=c(1), color="black", size=1) +
        scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0)) +
        scale_fill_gradientn(limits=c(0,outcome_threshold),
                             colours=rev(brewer.pal(11,"RdYlBu")),
                             values=c(seq(0,1,length.out=6),seq(1,outcome_threshold,length.out=6)[-1])/outcome_threshold,
                             breaks=c(0,1,2,outcome_threshold/2,outcome_threshold),
                             labels=c(0,1,2,outcome_threshold/2,outcome_threshold),
                             name="fold\nchange") +     
        xlab("average time of the direct switch, days") +
        ylab("average time of the inverse switch, days") +
        coord_equal() +
        theme_bw(base_size=base_sz,base_family='Times') +
#         guides(fill=FALSE) +
        annotate("point",x=28,y=60,shape=23,fill=cl_pt,size=2) +
        annotate("rect",xmin=27-w,xmax=27+w,ymin=105-h,ymax=105+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=27,y=105,label="4",size=3,vjust=vjst,hjust=hjst) +
        annotate("rect",xmin=62.5-w,xmax=62.5+w,ymin=105-h,ymax=105+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=62.5,y=105,label="2",size=3,vjust=vjst,hjust=hjst) +
        annotate("rect",xmin=105-w,xmax=105+w,ymin=67-h,ymax=67+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=105,y=67,label="1/2",size=2.5,vjust=vjst,hjust=hjst) +
        annotate("rect",xmin=105-w,xmax=105+w,ymin=45-h,ymax=45+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=105,y=45,label="1/4",size=2.5,vjust=vjst,hjust=hjst) +
        theme(legend.title=element_text(size=10,vjust=10),legend.text=element_text(size=8),
              axis.title=element_text(size=11)) -> p1
p1

ggsave(plot=p1,width=cs[1],height=cs[2],filename="../figures/draft/Fig5A.pdf",useDingbats=FALSE)<jupyter_output><empty_output><jupyter_text>Two other supplementary figures for $c=2\%$ and $c=6\%$<jupyter_code>cs = c(3.8,3)
options(repr.plot.width=cs[1],repr.plot.height=cs[2])

outcome_threshold = 10
nbin=20; w=3; h=2.2
cl_rc = brewer.pal(8,"Greys")[2]
cl_pt = "green" #brewer.pal(8,"Greens")[3]
vjst = .4
hjst = .5
df %>% filter(index==2) %>%
    mutate(outcome_cut=if_else(outcome>=outcome_threshold,outcome_threshold,outcome)) %>%
    ggplot(aes(x=T,y=Tbar,fill=outcome_cut)) + 
        geom_raster(interpolate=T) +
        geom_abline(slope=1,size=.2,linetype="dashed") + 
        stat_contour(aes(z=outcome),breaks=c(.2,.5,2,5), color=brewer.pal(8,"Greys")[1], size=.9) +
        stat_contour(aes(z=outcome),breaks=c(1), color="black", size=.9) +
        scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0)) +
        scale_fill_gradientn(limits=c(0,outcome_threshold),
                             colours=rev(brewer.pal(11,"RdYlBu")),
                             values=c(seq(0,1,length.out=6),seq(1,outcome_threshold,length.out=6)[-1])/outcome_threshold,
                             breaks=c(1e-4,1,2,outcome_threshold/2,outcome_threshold),
                             labels=c(0,1,2,outcome_threshold/2,outcome_threshold),
                             name="fold\nchange") + 
        xlab("average time of the direct switch, days") +
        ylab("average time of the inverse switch, days") +
        coord_equal() +
        theme_bw(base_size=base_sz,base_family='Times') +
        annotate("rect",xmin=105-w,xmax=105+w,ymin=95.5-h,ymax=95.5+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=105,y=95.5,label="4",size=3,vjust=vjst,hjust=hjst) +
        annotate("rect",xmin=105-w,xmax=105+w,ymin=61.5-h,ymax=61.5+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=105,y=61.5,label="2",size=3,vjust=vjst,hjust=hjst) +
        annotate("rect",xmin=105-w,xmax=105+w,ymin=37.5-h,ymax=37.5+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=105,y=37.5,label="1/2",size=2.5,vjust=vjst,hjust=hjst) +
        annotate("rect",xmin=105-w,xmax=105+w,ymin=28.5-h,ymax=28.5+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=105,y=28.5,label="1/4",size=2.5,vjust=vjst,hjust=hjst) +
        theme(legend.title=element_text(size=10,vjust=10),legend.text=element_text(size=8),
              axis.title=element_text(size=10)) -> p2
p2

ggsave(plot=p2,width=cs[1],height=cs[2],filename="../figures/draft/Fig5B.pdf",useDingbats=FALSE)
outcome_threshold = 10
nbin=20; w=3; h=2.2
cl_rc = brewer.pal(8,"Greys")[2]
cl_pt = "green" #brewer.pal(8,"Greens")[3]
vjst = .4
hjst = .5
df %>% filter(index==3) %>%
    mutate(outcome_cut=if_else(outcome>=outcome_threshold,outcome_threshold,outcome)) %>%
    ggplot(aes(x=T,y=Tbar,fill=outcome_cut)) + 
        geom_raster(interpolate=T) +
        geom_abline(slope=1,size=.2,linetype="dashed") + 
        stat_contour(aes(z=outcome),breaks=c(.2,.5,2,5), color=brewer.pal(8,"Greys")[1], size=.9) +
        stat_contour(aes(z=outcome),breaks=c(1), color="black", size=.9) +
        scale_fill_gradientn(limits=c(0,outcome_threshold),
                             colours=rev(brewer.pal(11,"RdYlBu")),
                             values=c(seq(0,1,length.out=6),seq(1,outcome_threshold,length.out=6)[-1])/outcome_threshold,
                             breaks=c(1e-4,1,2,outcome_threshold/2,outcome_threshold),
                             labels=c(0,1,2,outcome_threshold/2,outcome_threshold),
                             name="fold\nchange") + 
        scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0)) +
        theme_bw(base_size=base_sz,base_family='Times') +
        xlab("average time of the direct switch, days") +
        ylab("average time of the inverse switch, days") +
#         annotate("rect",xmin=27-w,xmax=27+w,ymin=105-h,ymax=105+h,fill=cl_rc,color="white",size=.6) +
#         annotate("text",x=27,y=105,label="4",size=3,vjust=vjst,hjust=hjst) +
        annotate("rect",xmin=7.5-w,xmax=7.5+w,ymin=105-h,ymax=105+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=7.5,y=105,label="2",size=3,vjust=vjst,hjust=hjst) +
        annotate("rect",xmin=45-w,xmax=45+w,ymin=105-h,ymax=105+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=45,y=105,label="1/2",size=2.5,vjust=vjst,hjust=hjst) +
        annotate("rect",xmin=101.5-w,xmax=101.5+w,ymin=105-h,ymax=105+h,fill=cl_rc,color="white",size=.6) +
        annotate("text",x=101.5,y=105,label="1/4",size=2.5,vjust=vjst,hjust=hjst) +
        coord_equal() +
        theme(legend.title=element_text(size=10,vjust=10),legend.text=element_text(size=8),
              axis.title=element_text(size=10)) -> p3
p3

ggsave(plot=p3,width=cs[1],height=cs[2],filename="../figures/draft/Fig5C.pdf",useDingbats=FALSE)<jupyter_output><empty_output><jupyter_text>Varying other parameters Cost of resistance<jupyter_code>"../results/final-sensitivity_c.csv" %>% 
    read.csv(colClasses='numeric') %>%
    select(-matches("err|sw_end|univ_point")) %>%
    arrange(c) -> df_

cs = c(3,2.75)
options(repr.plot.width=cs[1],repr.plot.height=cs[2])

outcome_bs = filter(df_,c==.04)$outcome

pC = df_ %>%
    gather(Variable,Value,-outcome) %>% mutate(Value = Value*100) %>%
    ggplot(aes(x=Value,y=outcome)) +
        geom_hline(yintercept=1,size=.25,color="black",linetype="dashed") +
        geom_hline(yintercept=0,size=.25,color="black",linetype="solid") +
        geom_line(size=.6) +
        coord_cartesian(expand=0.01,ylim=c(0,60)) +
        ylab("fold change in tumor size") + xlab("cost of resistance, %") +
        theme_bw(base_size=base_sz,base_family='Times') + 
        theme(
            panel.grid.major = element_line(colour="grey",size=.2), 
            panel.grid.minor = element_blank(),
            plot.margin=unit(c(.5,.5,.75,.5),"lines")) + 
        annotate("point",x=4,y=outcome_bs,shape=23,fill=cl_pt,size=2)

pC<jupyter_output><empty_output><jupyter_text>Sensitivity on $\alpha$<jupyter_code>"../results/final-sensitivity_alpha.csv" %>% 
    read.csv(colClasses='numeric') %>%
    select(-matches("err|sw_end|univ_point")) %>%
    arrange(alpha) -> df_

cs = c(3,2.75)
options(repr.plot.width=cs[1],repr.plot.height=cs[2])

pAlpha = df_ %>% gather(Variable,Value,-outcome) %>%
    ggplot(aes(x=Value,y=outcome)) + 
        geom_hline(yintercept=1,size=.25,color="black",linetype="dashed") +
        geom_hline(yintercept=0,size=.25,color="black",linetype="solid") +
        geom_line(size=.6) +
        coord_cartesian(expand=0.01,ylim=c(0,8)) +
        ylab("fold change in tumor size") + xlab(expression(alpha))  +
        theme_bw(base_size=base_sz,base_family='Times') + 
        theme(
            axis.title.y = element_text(colour=NA),
            panel.grid.major = element_line(colour="grey",size=.2), 
            panel.grid.minor = element_blank(),
            plot.margin=unit(c(.5,.5,.75,.9),"lines")) + 
        annotate("point",x=.3,y=outcome_bs,shape=23,fill=cl_pt,size=2)
pAlpha
"../results/final-sensitivity_L.csv" %>% 
    read.csv(colClasses='numeric') %>%
    select(-matches("err|sw_end|univ_point")) %>%
    arrange(L) -> df_

cs = c(3,2.75)
options(repr.plot.width=cs[1],repr.plot.height=cs[2])

pL = df_ %>% gather(Variable,Value,-outcome) %>%
    ggplot(aes(x=Value,y=outcome)) + 
        geom_hline(yintercept=1,size=.25,color="black",linetype="dashed") +
        geom_hline(yintercept=0,size=.25,color="black",linetype="solid") +
        geom_line(size=.6) +
        coord_cartesian(expand=0.01,ylim=c(0,3)) +
        ylab("fold change in tumor size") + xlab(expression(italic(L))) +
        theme_bw(base_size=base_sz,base_family='Times') + 
        theme(
            axis.title.y = element_text(vjust=4),
            panel.grid.major = element_line(colour="grey",size=.2), 
            panel.grid.minor = element_blank(),
            plot.margin=unit(c(.5,.5,.75,.75),"lines")) + 
        annotate("point",x=.2,y=outcome_bs,shape=23,fill=cl_pt,size=2)

pL
"../results/final-sensitivity_theta.csv" %>% 
    read.csv(colClasses='numeric') %>%
    select(-matches("err|sw_end|univ_point")) %>%
    arrange(theta) -> df_

cs = c(3,2.75)
options(repr.plot.width=cs[1],repr.plot.height=cs[2])

pTheta = df_ %>% gather(Variable,Value,-outcome) %>%
    ggplot(aes(x=Value,y=outcome)) + 
        geom_hline(yintercept=1,size=.25,color="black",linetype="dashed") +
        geom_hline(yintercept=0,size=.25,color="black",linetype="solid") +
        geom_line(size=.6) +
        coord_cartesian(expand=0.01,ylim=c(0,8),xlim=c(.25,.55)) +
        ylab("fold change in tumor size") + xlab(expression(theta)) +
        theme_bw(base_size=base_sz,base_family='Times') + 
        theme(
            axis.title.y = element_text(colour=NA),
            panel.grid.major = element_line(colour="grey",size=.2), 
            panel.grid.minor = element_blank(),
            plot.margin=unit(c(.5,.5,.75,.5),"lines")) + 
        annotate("point",x=.45,y=outcome_bs,shape=23,fill=cl_pt,size=2)

pTheta
"../results/final-sensitivity_kappa.csv" %>% 
    read.csv(colClasses='numeric') %>%
    select(-matches("err|sw_end|univ_point")) %>%
    arrange(kappa) -> df_

cs = c(3,2.75)
options(repr.plot.width=cs[1],repr.plot.height=cs[2])

pKappa = df_ %>% gather(Variable,Value,-outcome) %>%
    ggplot(aes(x=Value,y=outcome)) +
        geom_hline(yintercept=1,size=.25,color="black",linetype="dashed") +
        geom_hline(yintercept=0,size=.25,color="black",linetype="solid") +
        geom_line(size=.6) +
        coord_cartesian(expand=0.01,ylim=c(0,16),xlim=c(0,100)) +
        ylab("fold change in tumor size") + xlab(expression(kappa)) +
        theme_bw(base_size=base_sz,base_family='Times') + 
        theme(
            axis.title.y = element_text(colour=NA),
            panel.grid.major = element_line(colour="grey",size=.2), 
            panel.grid.minor = element_blank(),
            plot.margin=unit(c(.5,.5,.75,.5),"lines")) + 
        annotate("point",x=40,y=outcome_bs,shape=23,fill=cl_pt,size=2)

pKappa
cs = c(9.5,6.5)
options(repr.plot.width=cs[1],repr.plot.height=cs[2]) 
pFinal = grid.arrange(ggplotGrob(pC), ggplotGrob(pAlpha), ggplotGrob(pTheta), 
             ggplotGrob(pL), ggplotGrob(pKappa), nrow=2)

ggsave(plot=pFinal,height=cs[2],width=cs[1],dpi=200,filename="../figures/FigS10.pdf",useDingbats=FALSE)
print(paste("Number of recorded plots: ",length(plot_point_sizes),sep=""))
for (k in 1:length(plot_point_sizes)) {
    filename = paste("../figures/draft/",nm[k],".pdf",sep="")
    print(paste("(",k,") ",filename,sep=""))
    cs = plot_point_sizes[[k]]
    pdf(filename,width=cs[1],height=cs[2]);
    replayPlot(eval(as.name(nm[k])));
    dev.off() }<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19953645174884113, 'fraction_numerical': 0.04445849136114623, 'mean_word_length': 4.087173990710968, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 81, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '542397', 'license': 'MIT', 'n_tokens_mistral': 6133, 'n_tokens_neox': 5692, 'n_words': 462, 'path': 'scripts/.ipynb_checkpoints/C. Sensitivity analysis-checkpoint.ipynb', 'repo_name': 'aakhmetz/AkhmKim2019Scripts', 'size': '14238'}