library("reticulate")
use_python("/Users/ryan/anaconda3/bin/python")
esig=import("esig")

load(file="Example Data.rda")

#Note that your answers when running this code may differ slightly based on the randomness inherent
#to the k-fold cross validation procedure resulting in slightly different values for bandwidth and
#number of predictors.

X1=as.matrix(ExampleData[,1:100])
X2=as.matrix(ExampleData[,101:200])
X3=as.matrix(ExampleData[,201:300])
X4=as.matrix(ExampleData[,301:400])
Ymean=as.vector(ExampleData[,401])
Y=as.vector(ExampleData[,402])

##### Linear mFPCA #####

X=cbind(X1,X2,X3,X4)
k=K.Select.Lin(X,Y,PCA=F)

Intervals=matrix(0,nrow=nrow(X),ncol=2)
for(i in 1:nrow(X)){
  Interval=lin.mFPCA(X[-i,],Y[-i],X[i,],k=k,PCA=F,alpha=0.05)
  Intervals[i,1]=Interval$LB
  Intervals[i,2]=Interval$UB
}

LinearResults=cbind(Intervals,Ymean,Y)

TheoreticalCoverage=rep(0,400)
TheoreticalCoverage=pnorm((LinearResults[,2]-LinearResults[,3]),0,1)-pnorm((LinearResults[,1]-LinearResults[,3]),0,1)
mean(TheoreticalCoverage,na.rm=T)
#0.9562757
sd(TheoreticalCoverage)
#0.06252944

EmpiricalCoverage=mean((LinearResults[,1]<=LinearResults[,4])&(LinearResults[,4]<=LinearResults[,2]))
EmpiricalCoverage
#0.95

Length=rep(0,400)
Length=LinearResults[,2]-LinearResults[,1]
median(Length,na.rm=T)
#5.242564
mad(Length)
#0.036568

##### Local Linear mFPCA #####

X=cbind(X1,X2,X3,X4)
k=k.Select.ll(X,Y,PCA=F)
h=h.Select.ll(X,Y,K=k,PCA=F)

Intervals=matrix(0,nrow=nrow(X),ncol=2)
for(i in 1:nrow(X)){
  Interval=ll.mFPCA(X[-i,],X[i,],Y[-i],h=h,k=k,PCA=F,alpha=0.05)
  Intervals[i,1]=Interval$LB
  Intervals[i,2]=Interval$UB
}

LocalLinearResults=cbind(Intervals,Ymean,Y)

TheoreticalCoverage=rep(0,400)
TheoreticalCoverage=pnorm((LocalLinearResults[,2]-LocalLinearResults[,3]),0,1)-pnorm((LocalLinearResults[,1]-LocalLinearResults[,3]),0,1)
mean(TheoreticalCoverage,na.rm=T)
#0.9578936
sd(TheoreticalCoverage)
#0.03707635

EmpiricalCoverage=mean((LocalLinearResults[,1]<=LocalLinearResults[,4])&(LocalLinearResults[,4]<=LocalLinearResults[,2]))
EmpiricalCoverage
#0.95

Length=rep(0,400)
Length=LocalLinearResults[,2]-LocalLinearResults[,1]
median(Length,na.rm=T)
#4.747463
mad(Length)
#0.09824281

##### Partial Linear Signature #####

t=seq(0,1,length.out=100)
Z1=apply(X1,1,mean)
Z2=apply(X2,1,mean)
Z3=apply(X3,1,mean)
Z4=apply(X4,1,mean)
Z=cbind(Z1,Z2,Z3,Z4)

h=h.SelectPL(Y,Z)

SignatureX=matrix(0,n,31)
for(i in 1:n){
  data=cbind(t,X1[i,],X2[i,],X3[i,],X4[i,])
  SignatureX[i,]=esig$tosig$stream2sig(data,as.integer(3))
}
SignatureX=SignatureX[,8:31]

Intervals=matrix(0,nrow=nrow(X),ncol=2)
for(i in 1:nrow(X)){
  Interval=mPL.Sig(SignatureX[-i,],Z[-i,],SignatureX[i,],Z[i,],Y[-i],h=h,alpha=0.05)
  Intervals[i,1]=Interval$LB
  Intervals[i,2]=Interval$UB
}

PartialLinearResults=cbind(Intervals,Ymean,Y)

TheoreticalCoverage=rep(0,400)
TheoreticalCoverage=pnorm((PartialLinearResults[,2]-PartialLinearResults[,3]),0,1)-pnorm((PartialLinearResults[,1]-PartialLinearResults[,3]),0,1)
mean(TheoreticalCoverage,na.rm=T)
#0.9477502
sd(TheoreticalCoverage)
#0.07481955

EmpiricalCoverage=mean((PartialLinearResults[,1]<=PartialLinearResults[,4])&(PartialLinearResults[,4]<=PartialLinearResults[,2]))
EmpiricalCoverage
#0.95

Length=rep(0,400)
Length=PartialLinearResults[,2]-PartialLinearResults[,1]
median(Length,na.rm=T)
#5.308081
mad(Length)
#0.1477007