rm(list=ls())
library(Hmisc)

setwd('~/research/gTiger/Article_citsci/open_code')
source('settings/parameters.r')
source('settings/paths.r')
source('settings/functions.r')

setwd(path_to_data)

reporting_propensity_sampled_stanfit = readRDS('reporting_propensity_sampled_stanfit.Rds')

rp = reporting_propensity_sampled_stanfit$rp
bs = reporting_propensity_sampled_stanfit$bs
bsims = reporting_propensity_sampled_stanfit$bsims
user_intercepts = reporting_propensity_sampled_stanfit$user_intercepts

participation_days_vec = rp$participation_days

mean_participation_days = mean(rp$participation_days)
sd_participation_days = sd(rp$participation_days)

describe(participation_days_vec)

ndays = max(participation_days_vec)

propensity_func = function(participation_days, intercept, coefficients=bs) {
	participation_days = (participation_days - mean_participation_days)/(2*sd_participation_days)
	intercept + coefficients[1]*participation_days + coefficients[2]*(participation_days^2) + coefficients[3]*(participation_days^3)

}

propensity_func_prob = function(participation_days, intercept, coefficients=bs) {
	participation_days = (participation_days - mean_participation_days)/(2*sd_participation_days)
	odds = exp(intercept + coefficients[1]*participation_days + coefficients[2]*(participation_days^2) + coefficients[3]*(participation_days^3))
	prob = odds/(1+odds)
	return(prob)
}


make_propensity_curve_plot = function(){
	par(mar=c(4.5, 4.5, 0,0))
	curve(propensity_func_prob(x, bs[1], bs[2:4]) , from=min(participation_days_vec), to=max(participation_days_vec), xlab="participation time (days)", ylab="probability of sending report", type="n", ylim=c(0, .13), cex.axis=1.2, cex.lab=1.8)
for(i in 1:nrow(bsims)){
	curve(propensity_func_prob(x, bsims[i, 1], bsims[i,2:4]), from=min(participation_days_vec), to=max(participation_days_vec), add=TRUE, col="grey")
}
curve(propensity_func_prob(x, bs[1], bs[2:4]), from=min(participation_days_vec), to=max(participation_days_vec), add=TRUE)
}

pdf(paste0(path_to_figures, 'propensity_curve_participation_time.pdf'))
make_propensity_curve_plot()
dev.off()


png(paste0(path_to_figures, 'propensity_curve_participation_time.png'))
make_propensity_curve_plot()
dev.off()


pdf(paste0(path_to_figures, 'propensity_curve_user_intercepts.pdf'))
d = density(user_intercepts)
plot(d, main="", xlab="participant contribution to log odds of report")
polygon(d, col="#756bb1")
dev.off()

pdf(paste0(path_to_figures, 'propensity_curve_participation_time_sampled_user_means.pdf'))
curve(propensity_func_prob(x, user_intercepts[1], bs[2:4]) , from=min(participation_days_vec), to=max(participation_days_vec), xlab="participation time (days)", ylab="probability of sending report", ylim=c(0, .6), type="n")
for(i in 1:length(user_intercepts)){
	curve(propensity_func_prob(x, user_intercepts[i], bs[2:4]), from=min(participation_days_vec), to=max(participation_days_vec), add=TRUE, col="#7F0000")
}
dev.off()


pdf(paste0(path_to_figures, 'propensity_curve_log_odds_participation_time_sampled_user_means.pdf'))
curve(propensity_func(x, user_intercepts[1], bs[2:4]) , from=min(participation_days_vec), to=max(participation_days_vec), xlab="participation time (days)", ylab="log odds of sending report", type="n", ylim=c(-8, 0))
for(i in 1:length(user_intercepts)){
	curve(propensity_func(x, user_intercepts[i], bs[2:4]), from=min(participation_days_vec), to=max(participation_days_vec), add=TRUE, col="#7F0000")
}
dev.off()

library(rgl)

propensity_func_prob_3d = function(x, y) {
	y = (y - mean_participation_days)/(2*sd_participation_days)
	odds = exp(x + bs[2]*y + bs[3]*(y^2) + bs[4]*(y^3))
	prob = odds/(1+odds)
	return(prob)
}


d = density(user_intercepts)

densities = d$y
cp = floor(seq(1, length(densities), len=16))
dc = sapply(cp[1:length(cp)-1], function(i) sum(densities[i:i+1]))

myColorRamp <- function(colors, values) {
    v <- (values - min(values))/diff(range(values))
    x <- colorRamp(colors)(v)
    rgb(x[,1], x[,2], x[,3], maxColorValue = 255)
}

cols <- myColorRamp(c("blue", "red"), dc) 

r3dDefaults$windowRect <- c(0, 100, 1400, 1400) 


curve_3d_smaller_labels <- function(propensity_func_prob_3d, x_range=c(mean(user_intercepts), max(user_intercepts)), y_range=c(min(participation_days_vec), max(participation_days_vec)), col=cols){ 
	if (!require(rgl) ) {stop("load rgl")}
	xvec <- seq(x_range[1], x_range[2], len=15)
	yvec <- seq(y_range[1], y_range[2], len=15)
	fz <- outer(xvec, yvec, FUN=propensity_func_prob_3d)
	open3d()
	par3d(cex=1)
	persp3d( xvec, yvec, fz, col=col, xlab="", ylab="", zlab="", axes=F)
	axes3d(edge = 'bbox', box=F, expand = 1.03)
	mtext3d("participant motivation (random intercept)", edge=c('x-+'), line=3)
	mtext3d("probability", edge=c('z+-'), line=3)
	mtext3d("participation time (days)", edge=c('y+-'), line=3)
 }

curve_3d_smaller_labels(propensity_func_prob_3d)

writeWebGL(dir=file.path(paste0(path_to_figures, "webGL_propensity3D")), width=1000)


curve_3d <- function(propensity_func_prob_3d, x_range=c(mean(user_intercepts), max(user_intercepts)), y_range=c(min(participation_days_vec), max(participation_days_vec)), col=cols){ 
	if (!require(rgl) ) {stop("load rgl")}
	xvec <- seq(x_range[1], x_range[2], len=15)
	yvec <- seq(y_range[1], y_range[2], len=15)
	fz <- outer(xvec, yvec, FUN=propensity_func_prob_3d)
	open3d()
	par3d(cex=1.4)
	persp3d( xvec, yvec, fz, col=col, xlab="", ylab="", zlab="", axes=F)
	axes3d(edge = 'bbox', box=F, expand = 1.03, cex.axis=2)
	mtext3d("random intercept", edge=c('x-+'), line=3, cex=2)
	mtext3d("probability", edge=c('z+-'), line=4, cex=2)
	mtext3d("participation time (days)", edge=c('y+-'), line=3, cex=2)
 }

curve_3d(propensity_func_prob_3d)

rgl.postscript(paste0(path_to_figures, "propensity3d.pdf"), fmt="pdf")

rgl.postscript(paste0(path_to_figures, "propensity3d.svg"), fmt="svg")

