# Greg Francis
# 11 January 2019
graphics.off()  # clear all graphics
rm(list=ls(all=TRUE))  # clear all variables

Months <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

# load results data file
RawData<-read.csv(file="Results/Results.csv",header=TRUE)

# Remove practice and demographic trials
raw1 <- subset(RawData, RawData$TrialType!="Practice" )
ResponseData <-subset(raw1, raw1$TrialType!="Demographics" )

# Code TimeOfDay
TimeOfDay<-c()
# Convert timestamp to filename that matches image angles
fileNames <-c()
for(trial in c(1:length(ResponseData$Timestamp))){
	thisTS = paste(ResponseData$Timestamp[trial], "", sep="")
	# parse timestamp string appropriately (matches python code in SaveImages.py)
	# split the date with space
	thisTime = unlist(strsplit(thisTS,' ' ) )
	
	# find the month (and number), day, year
	thisMonth = thisTime[2]
	thisMonthNumber = 1 + match(Months, thisMonth)
	thisDate = as.numeric(thisTime[3])
	thisYear = thisTime[4]
	
	clockItems = unlist(strsplit(thisTime[5], ":") )
	hour = as.numeric(clockItems[1])

	minute = as.numeric(clockItems[2])
	second = as.numeric(clockItems[3])
	
	TimeOfDay<-c(TimeOfDay, hour+((minute+second/60)/60))
	
	if (thisDate <10){
		txtDate = paste("0",thisDate, sep="")
	} else{
		txtDate = thisDate
		}
	if (hour<10){
		txthour = paste("0",hour, sep="")
		}else {
		txthour = hour
		}
	if (minute < 10){
		txtminute=paste("0",minute, sep="")
		}	else{
		txtminute=minute	
		}
	if (second<10){
		txtsecond = paste("0",second, sep="")
		}	else{
		txtsecond =second
		}
		locale = "WestLafayette"

	fileName <- paste(locale,thisMonth,"-",txtDate,"-",thisYear,"h",txthour,"m",txtminute,"s",txtsecond,".jpg", sep="")
	fileNames<-c(fileNames, fileName)	
}
ResponseData$FileNames = fileNames
ResponseData$TimeOfDay = TimeOfDay

################
OrientationJudgments <-subset(ResponseData, ResponseData$TrialType=="Experiment")

# load altitude measures
ImageAltitudes<-read.csv(file="Results/Altitudes.csv",header=TRUE)

# Find altitude for each trial
altitudes<-c()
for(name in ResponseData$FileNames){
	responses <- subset(ImageAltitudes, ImageAltitudes$Image_filename==name)  # Should only be one match
	thisMean <- mean(responses$AltitudeDecimalDegrees)
	altitudes <-c(altitudes, thisMean)
}

OrientationJudgments$Altitude <-altitudes
#############################

# load TrueMoon angles (orientation of marker relative to straight down)
TrueMoonAngles<-read.csv(file="Results/TrueMoonOrientation.csv",header=TRUE)

# Average true angle for each trial (0 is straight down, postive runs counterclockwise to 360)
angles<-c()
MoonMarker<-c()
for(name in ResponseData$FileNames){
	responses <- subset(TrueMoonAngles, TrueMoonAngles$Image_filename==name)
	thisMean <- mean(responses$Alpha)
	# if negative (clockwise), then convert to positive (counterclockwise)
	if(thisMean<0) {thisMean = 360 + thisMean}	
	angles<-c(angles, thisMean)
	MoonMarker<-c(MoonMarker, responses$FeatureIdentified[1])  
}
OrientationJudgments$TrueMoonAngle <-angles
OrientationJudgments$MoonMarker <- MoonMarker

#########################
# load PhaseMoon angles (orientation of marker relative to straight down)
PhaseMoonAngles<-read.csv(file="Results/PhaseMoonOrientations.csv",header=TRUE)

# Average true angle for each trial (0 is straight down, postive runs counterclockwise to 360)
angles<-c()
for(name in ResponseData$ImageOnPhone){
	responses <- subset(PhaseMoonAngles, PhaseMoonAngles$Image_filename==name)
	thisMean <- mean(responses$Alpha)
	if(thisMean<0) {thisMean = 360 + thisMean}
	angles<-c(angles, thisMean)
}
# Make orientation be between -180 and 180
OrientationJudgments$PhaseMoonAngle <-angles

####################
# Compute orientation of match moon (relative to default PhaseMoon image)
# From the experiment MatchOrientation runs positive->Clockwise, negative->Counterclockwise (need it to be the other way around)
# Also want to run from 0 to 360
matchAngles <- -OrientationJudgments$MatchOrientation_size  # flip direction (CW to CCW)
matchAngles[matchAngles<0] <- 360 + matchAngles[matchAngles<0] # Runs from 0 to 360 CCW
OrientationJudgments$MatchOrientation <- matchAngles

judgedAngle <- (OrientationJudgments$MatchOrientation + OrientationJudgments$PhaseMoonAngle)
judgedAngle[judgedAngle <0] <- 360 + judgedAngle[judgedAngle <0] # Runs from 0 to 360 CCW
judgedAngle <- judgedAngle%%360 

OrientationJudgments$judgedAngle <-judgedAngle
OrientationJudgments$Errors <-((OrientationJudgments$judgedAngle - OrientationJudgments$TrueMoonAngle +180)%%360 -180 )

OrientationJudgments$sqErrors <-(OrientationJudgments$Errors)^2

write.csv(OrientationJudgments, file = "OrientationJudgmentResults.csv")


# Uncomment the following two lines to explore data for just waning (Aristarchus) or waxing (MareCrisium) moons
#OJ2<- subset(OrientationJudgments, OrientationJudgments$MoonMarker==2)  # 1 for Aristarchus, 2 for MareCrisium
#OrientationJudgments <- OJ2

#######################
# Remove outliers
NoOutliers<- subset(OrientationJudgments, OrientationJudgments$Errors < 45)
NoOutliers<- subset(NoOutliers, NoOutliers$Errors > -45)

cat("Number of trials before outlier removal: ", length(OrientationJudgments$Errors), "\n")
cat("Number of trials after outlier removal: ", length(NoOutliers$Errors), "\n")

# Compute average abs error for each subject
meanErrors <-aggregate(NoOutliers$Errors, list(NoOutliers $Participant), mean) 
meanAbsErrors <-abs(aggregate(NoOutliers$Errors, list(NoOutliers $Participant), mean) )
meanAltitudes <- aggregate(NoOutliers$Altitude, list(NoOutliers $Participant), mean)
maxAltitudes <- aggregate(NoOutliers$Altitude, list(NoOutliers $Participant), max)
minAltitudes <- aggregate(NoOutliers$Altitude, list(NoOutliers $Participant), min)
rangeAltitudes <- maxAltitudes$x - minAltitudes$x

# Find time spent on experiment
startTime <-aggregate(NoOutliers $TimeOfDay, list(NoOutliers $Participant), min) 
stopTime <-aggregate(NoOutliers $TimeOfDay, list(NoOutliers $Participant), max) 


cat("Number of subjects in full data set = ", length(unique(OrientationJudgments$Participant)),"\n")

cat("Number of subjects after outlier removal = ", length(meanAbsErrors$x),"\n")

cat("Correlation: ", cor(meanAltitudes$x, meanAbsErrors$x), "\n")

print (cor.test(meanAltitudes$x, meanAbsErrors$x))

cat("Done with frequentist analyses. \n")



#########
# Bayesian analysis
MIdata<- data.frame(AbsError=meanAbsErrors$x, Altitude = meanAltitudes$x)

# load the rethinking library
library(rethinking)

# fit a linear model that predicts absolute error as a function of moon elevation
MIPositiveSlope <- ulam(
	alist( AbsError ~ dnorm(mu, sigma), 
	mu <- a0 + a1*Altitude,
	a0 ~ dnorm(10, 10),
	a1 ~ dnorm(0, 20),
	sigma  ~ dunif(0, 100)
	), data= MIdata, constraints=list(a1="lower=0,upper=100"), log_lik=TRUE, iter=10000, warmup=4000, chains=4, cores=4  ) 

MILinear <- ulam(
	alist( AbsError ~ dnorm(mu, sigma), 
	mu <- a0 + a1*Altitude,
	a0 ~ dnorm(10, 10),
	a1 ~ dnorm(0, 20),
	sigma  ~ dunif(0, 100)
	), data= MIdata, log_lik=TRUE, iter=10000, warmup=4000, chains=4, cores=4  ) 

MINull <- ulam(
	alist( AbsError ~ dnorm(mu, sigma), 
	mu <- a0,
	a0 ~ dnorm(10, 10),
	sigma  ~ dunif(0, 100)
	), data= MIdata, log_lik=TRUE, iter=10000, warmup=4000, chains=4, cores=4  ) 


print(compare(MILinear, MIPositiveSlope, MINull) )



# Look at posterior distribution

samples <- extract.samples(MILinear)

cat("Probability slope is positive = ",  length(samples$a1[samples$a1 >0])/length(samples$a1) )


cat("All done!\n")