### CHANGE ALL PARTS BEFORE RUNNING CODE
#<your directory>
#<file ref>

library(fdasrvf)

start.time <- Sys.time()

#### CHANGE HERE - Find the file path of the lp points file (peat only)

path1root<-"<your directory>"

path3LPpoints<-paste("LP_points_20_modified_<file ref>.txt")
path3npeatM<-paste("MeanSeries_nopeat_<file ref>.csv")

path4out<-paste("ConditionMetrics_times_<file ref>_vA4_b.csv")
path5out<-paste("ConditionMetrics_<file ref>_vA4_b.csv")

Z<-read.table(paste(path1root,path3LPpoints,sep=""),skip=1)

Timev<-read.table(paste(path1root,path3LPpoints,sep=""),nrows=1)
Time.vecS2<-Timev-min(Timev)+1

mnpeat<-as.matrix(read.csv(paste(path1root,path3npeatM,sep=""),sep=","))
mnpeat<-as.numeric(mnpeat)

out1<-smooth.spline(Time.vecS2,mnpeat,spar=1)
out.pred1<-predict(out1,seq(min(Time.vecS2),max(Time.vecS2),by=6))
smooth.mat1<-out.pred1$y

mnpVel<-diff(smooth.mat1,differences=1)


###### FROM HERE (TO COPY)

k<-dim(Z)[2]

ID.vec1<-Z[c(148484:296966),1]
Lat.vec <- Z[c(148484:296966),2]
Long.vec  <- Z[c(148484:296966),3]
Z<-Z[c(148484:296966),-c(1:7)]


#### 1. CUBIC SPLINES TO SPLIT OSCILLATIONS FROM TRENDS

TimeMax<-as.numeric(max(Time.vecS2)) # number of daily time points between start and end date.
p<-dim(Z)[1] # number of peatland locations.


smooth.mat<-c()
smooth.mat1<-c()
smooth.mat2 <- c()

velf4<-c()
vell4<-c()

osc.vec3f4<-c()
osc.vec3l4<-c()


TimeMax<-max(Time.vecS2)
timeseq<-seq(1,TimeMax,by=6)

### CHECK SINE TEMPLATE

peakshift<-as.numeric(221)

sinfunc<-0.002*sin(2*pi*((timeseq+peakshift)/365))
sinfuncf4<-sinfunc[1:243]
sinfuncl4<-sinfunc[(467-242):467]

qsinef4 <- f_to_srvf(sinfuncf4,time=timeseq[1:243])
qsinel4 <- f_to_srvf(sinfuncl4,time=timeseq[(467-242):467])

eucliddist1.pl.f4<-function(x){sum(abs(x-qsinef4))}
eucliddist1.pl.l4<-function(x){sum(abs(x-qsinel4))}


for(i in 1:p){
  id <-i
  Z2<-Z[id,]
  Time.vecT22<-Time.vecS2
  if(any(is.na(Z[id,]==TRUE))==TRUE){
    Z2<-as.numeric(Z[id,-which(is.na(Z[id,]))])
    Time.vecT22<-as.numeric(Time.vecS2[-which(is.na(Z[id,]))])
  }
  Z3<-Z2
  Time.vecT23<-Time.vecT22
  k<-which(Z2==(-Inf))
  if(length(k)>=1){
    Z3<-as.numeric(Z2[-k])
    Time.vecT23<-as.numeric(Time.vecT22[-k])
  }
  Z4<-Z3
  Time.vecT24<-Time.vecT23
  k<-which(Z3==(Inf))
  if(length(k)>=1){
    Z4<-as.numeric(Z3[-k])
    Time.vecT24<-as.numeric(Time.vecT23[-k])
  }
  
  if(TimeMax!=max(Time.vecT24)){
    velf4[i]<-"NA"
    vell4[i]<-"NA"
    osc.vec3f4[i]<-"NA"
    osc.vec3l4[i]<-"NA"
  }
  
  if(TimeMax==max(Time.vecT24)){
  out1<-smooth.spline(Time.vecT24,Z4,spar=0.7)
  out.pred1<-predict(out1,seq(min(Time.vecT24),max(Time.vecT24),by=6))
  out2<-smooth.spline(Time.vecT24,Z4,spar=1)
  out.pred2<-predict(out2,seq(min(Time.vecT24),max(Time.vecT24),by=6))
  smooth.mat2<-out.pred2$y
  smooth.mat<-out.pred1$y-out.pred2$y
  
  nsub<-length(smooth.mat2)-1
  d2<-diff(smooth.mat2,differences=1)-mnpVel
  
  velf4[i]<-mean(d2[1:243])
  vell4[i]<-mean(d2[(nsub-242):nsub])
  
  qfunc <- f_to_srvf(smooth.mat[1:243],time=timeseq[1:243])
  osc.vec3f4[i]<-eucliddist1.pl.f4(qfunc)
  
  qfunc <- f_to_srvf(smooth.mat[(467-242):467],time=timeseq[(467-242):467])
  osc.vec3l4[i]<-eucliddist1.pl.l4(qfunc)
  
  }

  cat("\r",c("time series",i, " of ", p ))
}

XT1<-cbind(ID.vec1,Lat.vec,Long.vec)
XT2<-cbind(velf4,vell4)
XT5<-cbind(osc.vec3f4,osc.vec3l4)

XT<-cbind(XT1,XT2,XT5)

end.time1 <- Sys.time()

write.csv(XT,file=path5out,row.names=F)

time.taken1<-end.time1 - start.time
write.csv(c(p,time.taken1),file=path4out,row.names=F)