#Palaeoverse 2024 how to do trace fossil analyses
library(fitdistrplus)
library(mclust)
library(methods)
library(plyr)
library(selectspm)
library(spatstat)
library(spatstat.utils)
library(VGAM)
library(XML)
library(rlist)

library(ecespa)

#input data in. svg file with lines on
#run dex filaments from Katie Delahooke: https://github.com/kmdelahooke/dex/blob/main/dex_filaments.R

dataSh<-dex_filaments(file.choose()) #


# Plot points for Ediacaran window creation. 
plot(dataSh[,c(1,3)],ylim=c(0,1500),xlim=c(0,1500))
points(dataSh[,c(2,4)],ylim=c(0,1500),xlim=c(0,1500))
win.tr.final<-clickpoly(add=TRUE)

#create psp
psp.sh<-psp(dataSh[,1],dataSh[,3],dataSh[,2],dataSh[,4],win.tr.final)

#create ppp midpoints
ppp.sh.len<-ppp(coords(midpoints.psp(psp.sh))[,1],coords(midpoints.psp(psp.sh))[,2],window=win.tr.final, marks=lengths_psp(psp.sh))
ppp.sh.ang<-ppp(coords(midpoints.psp(psp.sh))[,1],coords(midpoints.psp(psp.sh))[,2],window=win.tr.final, marks=angles.psp(psp.sh))

#nice plots #using automatation to determine the appropraite density. 
tr.den<-density(psp.sh,bw.diggle(midpoints.psp(psp.sh)))
plot(tr.den,main="")
plot(psp.sh,add=TRUE)

#experiment
plot(density(ppp.sh.ang,10))

#list of the Segment lengths (each trace). 
tr.sh.len<-(lengths_psp(psp.sh))


par(mfrow=c(1,2))
hist(marks(ppp.sh.len)*100,main="Trace length",xlab="Segment length (cm)",col="grey",xaxs="i",yaxs="i",font.lab=2)
rose(marks(ppp.sh.ang)*180/pi,breaks=36,col="grey",main="Trace angle")


#### midpoint and intersection pcf analsyes
plot(envelope(ppp.sh.ang,pcf,nsim=99, nrank=4),legend=FALSE,xaxs="i",yaxs="i",font.lab=2,main="pcf",ylim=c(0,3),xlab="Length (m)",ylab="Midpoint PCF")

#with smoothing
plot(envelope(ppp.sh.ang,pcf,nsim=99, nrank=4,funargs = c(stoyan=0.9)),legend=FALSE,xaxs="i",yaxs="i",font.lab=2,main="pcf stoyan=0.9",ylim=c(0,3),xlab="Length (m)",ylab="Midpoint PCF")

#comparisons with CSR
 
fit.sh.an2<-envelope(ppm(unmark(ppp.sh.ang), trend=~1,statistic="pcf"),Gest,savefuns=TRUE,nsim=999,rank=50))

LF.gof(fit.sh.an2)

##### mark correlations functions

env.ppp.sh.ang<-envelope(ppp.sh.ang,markcorr,nsim=99,nrank=4,savefuns=TRUE)
env.ppp.sh.len<-envelope(ppp.sh.len,markcorr,nsim=99,nrank=4,savefuns=TRUE)

plot(env.ppp.sh.ang,legend=FALSE,xaxs="i",yaxs="i",font.lab=2,main="MCF big angles",xlab="Length (m)",ylab="MCF")
plot(env.ppp.sh.len,legend=FALSE,xaxs="i",yaxs="i",font.lab=2,main="MCF big lengths",xlab="Length (m)",ylab="MCF")
 
LF.gof(env.ppp.sh.ang)#$p#,ginterval=c(0,50)
LF.gof(env.ppp.sh.len)#$p#,ginterval=c(0,50)



par(mfrow=c(2,1))
plot(env.ppp.sh.len,legend=FALSE,xaxs="i",yaxs="i",font.lab=2,main="MCF lengths",xlab="Length (cm)",ylab="MCF",ylim=c(0.5,1.5))
plot(markcorr(ppp.sh.len),lty=2,add=TRUE)

plot(env.ppp.sh.ang,legend=FALSE,xaxs="i",yaxs="i",font.lab=2,main="MCF angles",xlab="Length (cm)",ylab="MCF",ylim=c(0.5,1.5))
plot(markcorr(ppp.sh.ang),lty=2,add=TRUE)

#Supplementary Work homogeneity tests
 
#find LH*
homtest(unmark(ppp.sh.ang),nsim=999)
sum(Gest(ppp.sh.ang)$han-Gest(ppp.mp.tr.f.big.angles)$theo)

descdist(tr.sh.len)

#best fit levy distribution
Coef(vglm(tr.sh.len ~ 1, levy(location=-0.01)))

#testing the distriubtions.
shapiro.test(tr.sh.len[is.finite((tr.sh.len))])

#comparison with best fit gaussian model (all very bad)
AIC(vglm(tr.sh.len ~ 1, levy(location=-0.01))) - AIC(glm(tr.sh.len ~ 1,family="gaussian"))

#testing the distriubtions.
shapiro.test(tr.sh.len[is.finite((tr.s))])