# This file contains code to generate the figures in
# An Empirical Model of Slope Ratio Comparisons
# Justin Talbot, John Gerth, and Pat Hanrahan
# InfoVis 2012

# This code tries to write out to a figures subdirectory, so you may need to create that.

library(ggplot2)
library(reshape2)
library(proto)
library(grid)
library(plyr)
library(scales)

# Common formatting

common.formatting <- 
( opts(panel.grid.minor=theme_blank(), panel.grid.major=theme_blank())
+ opts(panel.background = theme_rect(fill='white', colour='gray80'))
+ opts(strip.background = theme_rect(fill='gray85', colour='gray85'))
+ opts(axis.ticks = theme_segment(colour = 'gray80', size = 0.25))
+ opts(legend.position="none")
+ opts(strip.text.x = theme_text(size = 8))
+ opts(strip.text.y = theme_text(size = 8,angle=-90))
+ opts(axis.title.x=theme_text(vjust=-0.4,family="Helvetica",size=8),axis.title.y=theme_text(vjust=0.15,angle=90,family="Helvetica",size=8))
+ opts(axis.text.x=theme_text(family="Helvetica",size=8,colour='black'),axis.text.y=theme_text(family="Helvetica",size=8,colour='black'))
)

### Section 3

# Create a custom plot type to show the angles in a scatter plot

angleGrob <- function(x, y, angle1, angle2, length=0.07, size=0.5, 
	colour="black", linetype=1, arrow=NULL){
	
   o <- 0.025
   segmentsGrob(
		x0= c(x-o, x-o/2), 
		y0= c(y-o, y-o),
		x1= c(x-o + length*cos(angle1), x -o/2 + length*cos(angle1)),
		y1= c(y-o + length*sin(angle1), y-o + length*tan(angle2)*cos(angle1)),
		default.units="native",
	        gp = gpar(col=colour, 
                         lwd=0.5, 
                          lty=linetype, 
                          lineend = "butt"), 
	        arrow = arrow
               )
}

GeomAngle <- proto(ggplot2:::Geom, {
  objname <- "angle"

  draw_groups <- function(., ...) .$draw(...)
  draw <- function(., data, scales, coordinates, na.rm = FALSE, ...) {    
    data <- remove_missing(data, na.rm, 
      c("x", "y", "t1", "t2", "size", "shape"), name = "geom_angle")
    if (empty(data)) return(zeroGrob())
    
    with(coord_transform(coordinates, data, scales), 
      ggname(.$my_name(), angleGrob(x, y, t1, t2))
    )
  }

  draw_legend <- function(., data, ...) {
  }

  icon <- function(.) {
    pos <- seq(0.1, 0.9, length=6)
    pointsGrob(x=pos, y=pos, pch=19, gp=gpar(col="black", cex=0.5), default.units="npc")
  }
  
  default_stat <- function(.) StatIdentity
  required_aes <- c("x", "y", "t1", "t2")
  default_aes <- function(.) aes(shape=16, colour="black", size=2, fill = NA, alpha = 1)
  
})

geom_angle <- function (mapping = NULL, data = NULL, stat = "identity", position = "identity", na.rm = FALSE, ...) {
	GeomAngle$new(mapping = mapping, data = data, stat = stat, position = position, na.rm = na.rm, ...)
}

# Our stimuli
p2 <- read.csv("results_study1.csv", header=TRUE)
p <- ddply(p2, .(theta1,theta2), nrow)

# Cleveland stimuli
cleveland.theta1 <- rep(c(55,35.5,19.7,10.1)/180*pi, each=11)
cleveland.theta2 <- atan(tan(cleveland.theta1)*seq(0.5,1,length=11))
cleveland <- as.data.frame(cbind(theta1=cleveland.theta1, theta2=cleveland.theta2))

postscript(file="figures/figure2.eps",
        paper="special",
        width=3.3,
        height=3.15,
        horizontal=FALSE,
        fonts="Helvetica",
        colormodel="rgb")

(ggplot(p, aes(x=((theta1+theta2)/2)/pi*180,y=tan(theta2)/tan(theta1))) 
+ geom_point(data=cleveland, size=1.0, color="steelblue", shape=1)
+ geom_angle(aes(t1 = theta1, t2=theta2, length=0.9))

+ opts(aspect.ratio=1)
+ opts(plot.margin = unit(c(0,0,0.25,0.4), "lines"))
+ scale_y_continuous(limits=c(-0.01,1.01),name=substitute(paste("Slope ratio (",italic(p[ij]),")")),breaks=c(0.11, 0.17, 0.25, 0.35, 0.48, 0.65, 0.87), labels=c("11%", "17%", "25%", "35%", "48%", "65%", "87%"))
+ scale_x_continuous(limits=c(0,90),name=substitute(paste("Mid-angle (",italic(theta[m]),")")), breaks=c(9,18,30,45,60,72,81), labels=c(expression(9*degree),expression(18*degree),expression(30*degree),expression(45*degree),expression(60*degree),expression(72*degree),expression(81*degree)))

+ common.formatting
)

dev.off()


### Section 6 (Experiment 1)

n <- read.csv("results_study1.csv", header=TRUE)
n$strategy <- relevel(n$strategy, ref="Height")

n$midangle <- ((n$theta1+n$theta2)/2)/pi*180

n$error <- n$response-n$actual
n$abs.error <- abs(n$error)

n$our.model.nojudgmenterror <- 
    ifelse(n$strategy=="Angle", n$theta2/n$theta1*100, n$actual)
n$our.model <- 
    ifelse(n$strategy=="Angle", n$theta2/n$theta1*100+(n$midangle*-0.40+23.8), n$actual + 7.45)

# Absolute error predicted by our model (to be on the same scale as Cleveland's model)
n$our.model.nojudgmenterror.abs.error <- abs(n$our.model.nojudgmenterror - n$actual)
n$our.model.abs.error <- abs(n$our.model - n$actual)

# Cleveland's model directly predicts absolute error
n$cleveland.model.abs.error <- (4.39 + -0.47*(n$actual-100) + -1.14*(n$theta1-n$theta2)/pi*180)

# Residuals
n$residuals.nojudgmenterror <- n$response-n$our.model.nojudgmenterror
n$residuals <- n$response-n$our.model

n$sr <- cut(n$actual, breaks=c(0,15,20,30,40,60,70,100), labels=c("11%", "17%", "25%", "35%", "48%", "65%", "87%"))
n$inbounds <- ifelse(as.integer(n$strategy=="Angle" | n$actual < 40 | n$midangle > 45), 1, 2)


postscript(file="figures/figure5.eps",
        paper="special",
        width=7,
        height=4,
        horizontal=FALSE,
        fonts="Helvetica",
	colormodel="rgb")

(ggplot(n, aes(x=midangle,abs.error)) 
+ facet_grid(strategy~sr)
+ stat_smooth(fill="gray90",method="loess",se=TRUE,lwd=0,alpha=1,span=0.6,colour="gray90")
+ geom_hline(yintercept=0, color="gray80", lwd=0.25)
+ geom_vline(xintercept=45, color="gray80", lwd=0.25)
+ geom_point(colour="black",size=0.6,position=position_jitter(w=0,h=1)) 
+ stat_smooth(lwd=0.35,method="loess",se=FALSE,colour="gray30", span=0.6)
+ geom_line(aes(x=midangle,y=our.model.nojudgmenterror.abs.error), colour="red", lwd=0.5)
+ scale_size(range=c(0.35,0.95))

+ opts(plot.margin = unit(c(0.25,0,1,1), "lines"))
+ scale_y_continuous(limits=c(-15,70), name="Absolute error (percentage points)", breaks=c(0,50), labels=c("0", "50"))
+ scale_x_continuous(limits=c(5,85), name=substitute(paste("Mid-angle (",italic(theta[m]),")")), breaks=c(10,45,80), labels=c(expression(10*degree),expression(45*degree),expression(80*degree)))

+ common.formatting
)

dev.off()


postscript(file="figures/figure6.eps",
        paper="special",
        width=2.5,
        height=3,
        horizontal=FALSE,
        fonts="Helvetica",
	colormodel="rgb")

(ggplot(n, aes(x=midangle,y=residuals.nojudgmenterror)) 
+ facet_grid(.~strategy)
+ stat_smooth(fill="gray80",colour="gray80",method="loess",se=TRUE,lwd=0,alpha=1)
+ geom_hline(yintercept=0, color="gray80", lwd=0.25)
+ geom_point(colour="black",size=0.6) 
+ stat_smooth(lwd=0.35,method="loess",se=FALSE,colour="gray30", span=0.6)

+ opts(plot.margin = unit(c(0,0.2,1,1), "lines"))
+ scale_y_continuous(limits=c(-30,70), name="Residuals (percentage points)", breaks=c(0,50), labels=c("0", "50"))
+ scale_x_continuous(limits=c(5,85), name=substitute(paste("Mid-angle (",italic(theta[m]),")")), breaks=c(10,45,80), labels=c(expression(10*degree),expression(45*degree),expression(80*degree)))

+ common.formatting
)

dev.off()


postscript(file="figures/figure7.eps",
        paper="special",
        width=7,
        height=4,
        horizontal=FALSE,
        fonts="Helvetica",
	colormodel="rgb")

(ggplot(n, aes(x=midangle,y=abs.error)) 
+ facet_grid(strategy~sr)
+ stat_smooth(fill="gray90",method="loess",se=TRUE,lwd=0,alpha=1,span=0.6,colour="gray90")
+ geom_hline(yintercept=0, color="gray80", lwd=0.25)
+ geom_vline(xintercept=45, color="gray80", lwd=0.25)
+ geom_point(colour="black",size=0.6,position=position_jitter(w=0,h=1)) 
+ stat_smooth(lwd=0.35,method="loess",se=FALSE,colour="gray30", span=0.6)
+ geom_line(aes(x=midangle,y=our.model.abs.error), colour="red", lwd=0.5)
+ geom_line(aes(x=midangle,y=cleveland.model.abs.error, size=inbounds), color="cornflowerblue")
+ scale_size(range=c(0.35,0.95))

+ opts(plot.margin = unit(c(0.25,0,1,1), "lines"))
+ opts(title = substitute(paste("Slope ratio (",italic(p[ij]),")")), plot.title=theme_text(size=8,vjust=1)) 
+ scale_y_continuous(limits=c(-15,70), name="Absolute error (percentage points)", breaks=c(0,50), labels=c("0", "50"))
+ scale_x_continuous(limits=c(5,85), name=substitute(paste("Mid-angle (",italic(theta[m]),")")), breaks=c(10,45,80), labels=c(expression(10*degree),expression(45*degree),expression(80*degree)))

+ common.formatting
)

dev.off()


# This figure does not appear in the paper. It tests for a learning effect in Experiment 1.
postscript(file="figures/ex1_learning.eps",
        paper="special",
        width=7,
        height=2.2,
        horizontal=FALSE,
        fonts="Helvetica",
	colormodel="rgb")

names(n)

(ggplot(n, aes(x=X %% 147,y=abs.error)) 
+ facet_grid(.~strategy)
+ stat_smooth(fill="gray60",method="loess",se=TRUE,lwd=0,alpha=1)
+ geom_hline(yintercept=0, color="black", lwd=0.25)
+ geom_point(colour="gray70",size=1.01) 
+ stat_smooth(lwd=0.6,method="loess",se=FALSE,colour="black", span=0.6)

+ scale_y_continuous(name="Absolute Error (percentage points)", breaks=c(0,50,100), labels=c("0", "50", "100"))
+ scale_x_continuous(name="Stimuli presentation order")

+ common.formatting
)

dev.off()


### Section 7 (Experiment 2)

n <- read.csv("results_study2.csv", header=TRUE)
n$strategy <- relevel(n$strategy, ref="Height")
n$condition <- relevel(n$condition, ref="Without Baseline")

n$midangle <- ((n$theta1+n$theta2)/2)/pi*180

n$error <- n$response-n$actual
n$abs.error <- abs(n$error)

n$our.model.nojudgmenterror <- 
    ifelse(n$strategy=="Angle", n$theta2/n$theta1*100, n$actual)
n$our.model <- 
    ifelse(n$strategy=="Angle", n$theta2/n$theta1*100+(n$midangle*-0.40+23.8), n$actual + 7.45)

# Absolute error predicted by our model (to be on the same scale as Cleveland's model)
n$our.model.nojudgmenterror.abs.error <- abs(n$our.model.nojudgmenterror - n$actual)
n$our.model.abs.error <- abs(n$our.model - n$actual)

# Cleveland's model directly predicts absolute error
n$cleveland.model.abs.error <- (4.39 + -0.47*(n$actual-100) + -1.14*(n$theta1-n$theta2)/pi*180)

# Residuals
n$residuals.nojudgmenterror <- n$response-n$our.model.nojudgmenterror
n$residuals <- n$response-n$our.model

n$sr <- cut(n$actual, breaks=c(0,15,20,30,40,60,70,1), labels=c("11%", "17%", "25%", "35%", "48%", "65%", "87%"))
n$inbounds <- ifelse(as.integer(n$strategy=="Angle" | n$actual < 40 | n$midangle > 45), 1, 2)


postscript(file="figures/figure8.eps",
        paper="special",
        width=7,
        height=4.5,
        horizontal=FALSE,
        fonts="Helvetica",
	colormodel="rgb")

(ggplot(n[n$strategy=="Angle",], aes(x=midangle,y=error)) 
+ facet_grid(condition~sr)
+ stat_smooth(fill="gray85",method="loess",se=TRUE,lwd=0,alpha=1,span=0.6,colour="gray85")
+ geom_hline(yintercept=0, color="gray80", lwd=0.25)
+ geom_vline(xintercept=45, color="gray80", lwd=0.25)
+ geom_point(colour="black",size=0.6) 
+ stat_smooth(lwd=0.35,method="loess",se=FALSE,colour="gray30", span=0.6)
+ scale_size(range=c(0.35,0.95))

+ opts(plot.margin = unit(c(0.25,0,1,1), "lines"))
+ opts(title = substitute(paste("Slope ratio (",italic(p[ij]),")")), plot.title=theme_text(size=8,vjust=1)) 
+ scale_y_continuous(limits=c(-30,70), name="Observed error (percentage points)", breaks=c(0,50), labels=c("0", "50"))
+ scale_x_continuous(limits=c(5,85), name=substitute(paste("Mid-angle (",italic(theta[m]),")")), breaks=c(10,45,80), labels=c(expression(10*degree),expression(45*degree),expression(80*degree)))

+ common.formatting
)

dev.off()


postscript(file="figures/figure9.eps",
        paper="special",
        width=2.5,
        height=3,
        horizontal=FALSE,
        fonts="Helvetica",
	colormodel="rgb")

(ggplot(n[n$strategy=="Angle",], aes(x=midangle,y=residuals.nojudgmenterror)) 
+ facet_grid(.~condition)
+ geom_hline(yintercept=0, color="gray80", lwd=0.25)
+ geom_point(colour="black",size=0.6) 
+ stat_smooth(fill="gray80",colour="gray80",method="loess",se=TRUE,lwd=0,alpha=1)
+ stat_smooth(lwd=0.35,method="loess",se=FALSE,colour="gray30", span=0.6)

+ opts(plot.margin = unit(c(0,0.2,1,1), "lines"))
+ scale_y_continuous(limits=c(-30,70), name="Residuals (percentage points)", breaks=c(0,50), labels=c("0", "50"))
+ scale_x_continuous(limits=c(5,85), name=substitute(paste("Mid-angle (",italic(theta[m]),")")), breaks=c(10,45,80), labels=c(expression(10*degree),expression(45*degree),expression(80*degree)))

+ common.formatting
)

dev.off()


# This figure does not appear in the paper. It tests for a learning effect in Experiment 2.
postscript(file="figures/ex2_learning.eps",
        paper="special",
        width=7,
        height=4,
        horizontal=FALSE,
        fonts="Helvetica",
	colormodel="rgb")

(ggplot(n, aes(x=X %% 147,y=abs.error)) 
+ facet_grid(strategy~condition)
+ stat_smooth(fill="gray60",method="loess",se=TRUE,lwd=0,alpha=1)
+ geom_hline(yintercept=0, color="black", lwd=0.25)
+ geom_point(colour="gray70",size=1.01) 
+ stat_smooth(lwd=0.6,method="loess",se=FALSE,colour="black", span=0.6)

+ scale_y_continuous(name="Absolute Error (percentage points)", breaks=c(0,50,100), labels=c("0", "50", "100"))
+ scale_x_continuous(name="Stimuli presentation order")

+ common.formatting
)

dev.off()

