rm(list=ls())
library(data.table)
library(rstan)
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)

rplist = readRDS('reporting_propensity_sampled_stanfit.Rds')

rp = rplist$rp
bs = rplist$bs

trs = readRDS('user_locations.Rds')
users = readRDS('user_reporting.Rds')

participation_days_vec = rp$participation_days
ndays = max(participation_days_vec)
mean_participation_days = mean(rp$participation_days)
sd_participation_days = sd(rp$participation_days)
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)
}

users$any_reports = users$n_reports > 0

report_prob_step1 = mean(users$any_reports)
trs$propensity = report_prob_step1 * propensity_func_prob(trs$participation_days, bs[1], bs[2:4])

saveRDS(trs, 'user_locations_propensities.Rds')
