library(RandomFields)
library(fBasics)
library(stabledist)
library(gtools)
library(Rfast)

#install.packages("stabledist")
#cat("\014")
#set.seed(999)
#set.seed(NULL)

#Dimention d=1
h=0.02   #mesh size
# Observations (past) points
T2 = 30           ### upper limit of time points
T1 = 0            ### lower limit of time points
RK = 5            # Prediction interval
#RK = 0.5            # Prediction interval
n=10                 # number of prediction points
step2 = 1


#Model=1 > Gaussian
#Model=2 > Cauchy
#Model=3 > Levy
#Model=4 > Student

#ModelPred=0 ##(0-intepolation, 1- extrapolation)


for (Model in 1:1){
  for (ModelPred in c(0,1)){
    
    
Model=1
ModelPred=0

#ModelPred=0 ##(0-intepolation, 1- extrapolation)


#Initialization
filename1<-switch (Model,  "GaussTrajA-h002.csv", "CauchyTrajA-h002.csv", "LevyTrajA-h002.csv","StudentTrajA-h002.csv")
filename2<-switch (Model,  "2108GradGauss-h002.csv", "2108GradCauhy-h002.csv", "2108GradLevy-h002.csv","2108GradStudent-h002.csv")
if (Model==1){
  F0<-function(x) return(pnorm(x,mean=-0.48,sd=1.11)) 
  p0<-function(x) return(dnorm(x,mean=-0.48,sd=1.11))
}

if (Model==2){
  F0<-function(x) return(pcauchy(x,scale=1.26)) 
  p0<-function(x) return(dcauchy(x,scale=1.26))
}

if (Model==3){
  mu=0
  cc<-0.836
  F0<-function(x) 
  {
    if (sum(x<0)>0) print("Alert")
    return(2-2*pnorm(sqrt(cc/(x))))
  }
  p0<-function(x) return((2*pi/cc)^(-0.5)*(x)^(-1.5)*exp(-0.5*cc/(x-mu)))
}

if (Model==4){
  F0<-function(x) return(pt(x/14,df=0.67)) 
  p0<-function(x) return(dt(x/14,df=0.67))
}

F0(c(-101.96,2))
p0(c(-101.96,2))


T3=T2+RK          # upper limit of predicted points

TA=seq(from=T1, to=T3-h,by=h)    #All time points
NTA=length(TA)                 #Number of all time points
ITA=1:NTA                      #All indeces of time points

T0=seq(from=T1, to=T2-h,by=h)    #All time points of observations
NT0=length(T0)                 #Number of observations
IT0=1:NT0                      #All indeces of observations



if (ModelPred==0){
  step1<-round((RK/h)/n)
  # Prediction points
  TP<-seq(from=T2, to=T2+RK-h,by=h*step1)    #All time points of Prediction points
  NTP<-length(TP)                 #Number of Prediction points
  ITP<-ITA[NT0+(0:(n-1))*step1+1]
  TA[ITP]
  # extrapolated points
  TE=TA[-append(IT0,ITP)]       #All time points of extrapolation
  NTE=length(TE)                #Number of extrapolated points
  ITE=ITA[-append(IT0,ITP)]     #indecies of extrapolated points
  
}else
{
  step1<-round((RK/h)/n)
  # Prediction points
  TP<-seq(from=T2, to=T2+5*h*(n-1),by=5*h)    #All time points of Prediction points
  #TP=seq(from=T2, to=T2+h*(n-1),by=h)    #All time points of Prediction points
  NTP=length(TP)                 #Number of Prediction points
  ITP<-ITA[NT0+(0:(NTP-1))*5+1]
  TA[ITP]
  # extrapolated points
  TE<-TA[-append(IT0,ITP)]       #All time points of extrapolation
  NTE<-length(TE)                #Number of extrapolated points
  ITE<-ITA[-append(IT0,ITP)]     #indecies of extrapolated points
}





TrajA=read.csv(filename1)
MaxN=max(TrajA[,1])
MaxT=max(TrajA[,2])
TrajA=TrajA[ITA+MaxN-NTA,3]
Traj=TrajA[IT0]

#Prediction:


#Gradient functional
hjarr<-function(ti) return(seq(from=1-min(ti,ITP),to=NT0-1-max(ti,ITP),by=step2))
nt2<-NT0+min(ITP)-max(ITP)-1
XPRarr<-matrix(rep(0,n*nt2),nrow=n)
for (i in 1:nt2)
{
  XPRarr[,i]<-Traj[ITP+i-min(ITP)]
}

#Functional

PhiF<-function(type,lambda,Xthat,Xt,gmm,iar)
{
  if (type==1)
  {
    return(mean(2*pmax(F0(Xthat),F0(Xt))-F0(Xthat)-F0(Xt)))
  }
  if (type==21)
  {
    Fxhat<-F0(Xthat)
    return(c(mean(abs(Fxhat-F0(Xt))),mean(gmm*(1/3-pmax(Fxhat,Fxhat[iar])+Fxhat^2))))
  }
  if (type==2)
  {
    Fxhat<-F0(Xthat)
    return(mean(abs(Fxhat-F0(Xt))+gmm*(1/3-pmax(Fxhat,Fxhat[iar])+Fxhat^2)))
  }
  if (type==3)
  {
    Fxhat<-F0(Xthat)
    return(mean(2*pmax(Fxhat,F0(Xt))-1+gmm*(1/3-pmax(Fxhat,Fxhat[iar])+Fxhat^2)))
  }
  if (type==4)
  {
    #Xthat<-((lambda0^2)%*%XPRti)[1,]
    #Xt<-Traj[ti+hia]
    Fxhat<-F0(Xthat)
    lenxa<-length(Fxhat)
    #ix<-combinations(n=lenxa, r=2, set = TRUE, repeats.allowed = FALSE)
    return(mean(abs(Fxhat-F0(Xt))+gmm*(1/3+Fxhat^2-Fxhat))-gmm*total.dist(Fxhat)/lenxa^2)
  }
}




#Gradient j
NablaPhiFj<-function(type,lambda,Xthat,Xt,XBhat,gmm,Yj,YBhat){
if (type==1)
{
  return(XBhat*(2*(Xt<Xthat)-1)*p0(Xthat))
}
if (type==2)
{
  return(XBhat*(2*(Xt<Xthat)-1+2*gmm*F0(Xthat)-gmm*(Yj<Xthat))*p0(Xthat))
}
if (type==3)
  {
    return(XBhat*(2*(Xt<Xthat)+2*gmm*F0(Xthat)-gmm*(Yj<Xthat))*p0(Xthat)-YBhat*gmm*(Yj>Xthat)*p0(Yj))
}
if (type==31)
  {
    return()
  }  
if (type==4)
{

  Fxhat<-F0(Xthat)
  lenxa<-length(Yj)
#  ix<-combinations(n=length(Fxhat), r=2, set = TRUE, repeats.allowed = FALSE)
#  return(mean(2*pmax(Fxhat,F0(Xt))-Fxhat-F0(Xt)+gmm*(1/3+Fxhat^2-Fxhat/length(Fxhat)-2*pmax(Fxhat[ix[,1]],Fxhat[ix[,2]])/length(Fxhat))))
  return(XBhat*(2*(Xt<Xthat)-1+2*gmm*F0(Xthat)-gmm*sum(Yj<=Xthat)/lenxa)*p0(Xthat)-gmm*YBhat%*%((Yj>Xthat)*p0(Yj))/lenxa)
}
}

#Gradient mean
NablaPhiF<-function(type,lambda,Xthat,Xt,XBhat,gmm,Yj)
{
  if (type==1)
  {
    return(rowMeans(XBhat*((2*(Xthat<Xt)-1)*p0(Xthat))))
  }
if (type==2)
{
  Fxhat<-F0(Xthat)
  return(mean(pmax(Fxhat,F0(Xt))+gmm*(pmax(Fxhat,sample(Fxhat,n=length(Fxhat)))-Fxhat^2)))
}
if (type==3)
{
  Fxhat<-F0(Xthat)
  
  return(mean(pmax(Fxhat,F0(Xt))+gmm*(pmax(Fxhat,sample(Fxhat,n=length(Fxhat)))-Fxhat^2)))
}
}



#unconditional 
LambdaMin<-function(type1,ti,nl,lmax,gm1,lambda0) {
# type1=4
#  gm1=5
#  ti=ITE[150]
#  lmax=300
#    nl=(10/(10+1:lmax))^0.7

    
  hia<-hjarr(ti)
  XPRti<-XPRarr[,hia+min(ITP)]
  lambda1<-matrix(rep(0,n*(lmax+1)),ncol = n)
#  Gradj1<-lambda1
#  Gradj2<-lambda1
  PhiL1<-rep(0,lmax)
  PhiMin=1000
  #iar=sample(1:length(hia),size = length(hia))
 
  iar<-(1:length(hia)+ sample((1*length(hia)%/%3):(2*length(hia)%/%3),size =  length(hia),replace = TRUE))%%length(hia)+1
 
  
  lambda1[1,]<-lambda0
  if (Model==3) {
  PhiMin<-PhiF(type=type1,Xthat = ((lambda0^2)%*%XPRti)[1,],Xt=Traj[ti+hia],gmm=gm1,iar=iar)
  }else
  {
    PhiMin<-PhiF(type=type1,Xthat = (lambda0%*%XPRti)[1,],Xt=Traj[ti+hia],gmm=gm1,iar=iar)
  }
  
  hiar<-hia[iar]
  for (j in 1:n)
  {
    lambdat<-runif(n)
    lambdat<-lambdat/sum(lambdat)
    if (Model==3) {
      Pht<-PhiF(type=type1,Xthat = ((lambdat^2)%*%XPRti)[1,],Xt=Traj[ti+hia],gmm=gm1,iar=iar)
    }
    else{
    Pht<-PhiF(type=type1,Xthat = (lambdat%*%XPRti)[1,],Xt=Traj[ti+hia],gmm=gm1,iar=iar)
    }
    if (Pht<PhiMin) {
      lambda1[1,]<-lambdat
      PhiMin<-Pht
    }
  }
  for (j in 1:n)
  {
    lambdat<-rep(0,n)
    lambdat[j]<-1
    if (Model==3) {
      Pht<-PhiF(type=type1,Xthat = ((lambdat^2)%*%XPRti)[1,],Xt=Traj[ti+hia],gmm=gm1,iar=iar)
    }
    else{
      Pht<-PhiF(type=type1,Xthat = (lambdat%*%XPRti)[1,],Xt=Traj[ti+hia],gmm=gm1,iar=iar)
    }
    if (Pht<PhiMin) {
      lambda1[1,]<-lambdat
      PhiMin<-Pht
    }
  }
  PhiL1[1]<-PhiMin
  lambdamin<-lambda1[1,]
  #hri<-sample(hia,size=lmax,replace = TRUE)
  

  
  
  
  iar2<-sample(1:length(hia),size = lmax,replace = TRUE)
  hri<-hia[iar2]
  #hri<-sample(hia,size=lmax,replace = TRUE)
  #hri2<-sample(hia,size=lmax,replace = TRUE)
  hri2<-hiar[iar2]
  XtA<-Traj[ti+hia]
  XBhatA<-XPRarr[,min(ITP)+hia]
  XBhatR<-XPRarr[,min(ITP)+hiar]
  
 GradQ<-0
  
if (GradQ==0){  
  
  for (l in 1:(lmax-1))
  {
    XBhatl<-XPRarr[,min(ITP)+hri[l]]
    YBhatl<-XPRarr[,min(ITP)+hri2[l]]
    Xhatl<-(lambda1[l,]%*%XBhatl)[1,]
    Yl<-(lambda1[l,]%*%YBhatl)[1,]
    Xtl<-Traj[ti+hri[l]]
    XhatA<-(lambda1[l,]%*%XBhatA)[1,]
    PhiL1[l]<-PhiF(type=type1,Xthat = XhatA,Xt=XtA,gmm=gm1,iar=iar)
    
    if (PhiL1[l]<PhiMin) {
      lambdamin<-lambda1[l,]
      PhiMin<-PhiL1[l]
    }
    
    
    
    if (type1==4){
      Gradj<-NablaPhiFj(type=type1,lambda=lambda1[l,],Xthat=Xhatl,Xt=Xtl,XBhat=XBhatl,gmm=gm1,Yj=(lambda1[l,]%*%XPRti)[1,],YBhat=XPRti)/length(hia)
   #   lenxa=length((lambda1[l,]%*%XPRti)[1,])
  #    XBhatl*(-gm1*sum((lambda1[l,]%*%XPRti)[1,]<=Xhatl)/lenxa)*p0(Xhatl)-gm1*XPRti%*%(((lambda1[l,]%*%XPRti)[1,]>Xhatl)*p0((lambda1[l,]%*%XPRti)[1,]))/lenxa
    }
    else  {
    Gradj<-NablaPhiFj(type=type1,lambda=lambda1[l,],Xthat=Xhatl,Xt=Xtl,XBhat=XBhatl,gmm=gm1,Yj=Yl,YBhat=YBhatl)/length(hia)
    }
    
    #    Gradj1[l,]<-NablaPhiFj(type=type1,lambda=lambda1[l,],Xthat=Xhatl,Xt=Xtl,XBhat=XBhatl,gmm=gm1,Yj=Yl)
#    Gradj2[l,]<-NablaPhiFj(type=1,lambda=lambda1[l,],Xthat=Xhatl,Xt=Xtl,XBhat=XBhatl,gmm=gm1,Yj=Yl)
   

    
#    lambda1[l+1,]<-lambda1[l,] - nl[l]*Gradj
#    Xhatl<-(lambda1[l+1,]%*%XBhatl)[1,]
#    Yl<-(lambda1[l+1,]%*%YBhatl)[1,]
#    XhatA<-(lambda1[l+1,]%*%XBhatA)[1,]
#    PhiL1[l+1]<-PhiF(type=type1,Xthat = XhatA,Xt=XtA,gmm=gm1,iar=iar)
#    if (PhiL1[l+1]>PhiL1[l]) {
#      lambda1[l+1,]<-lambda1[l,]
#    }
#    if (PhiL1[l]<PhiMin) {
#      lambdamin<-lambda1[l,]
#      PhiMin<-PhiL1[l]
#    }
    
    
     lambda1[l+1,]<-lambda1[l,] - nl[l]*Gradj
    
    
    
    #    XhatA<-(lambda1[l,]%*%XBhatA)[1,]
#      Grad<-NablaPhiF(type=type1,lambda=lambda1[l,],Xthat=XhatA,Xt=XtA,XBhat=XBhatA,gmm=gm1,Yj=Yl)
#      lambda1[l+1,]<-lambda1[l,] -0.001*Grad
      }
  
  #plot( PhiL1,type="l")
 # print(c(PhiF(type=21,Xthat =(lambdamin%*%XBhatA)[1,],Xt=XtA,gmm=gm1,iar=iar)))
#  lambdaN<-c(0,0,0,0,0,0,0,0,0,1)
#  print(c("last",PhiF(type=21,Xthat =(lambdaN%*%XBhatA)[1,],Xt=XtA,gmm=gm1,iar=iar)))
}else
 {
  PhiFL<-function(x) return(100*PhiF(type=type1,lambda=x^2,Xthat=((x^2)%*%XBhatA)[1,],XtA,gmm=gm1,iar=iar))
  ROU=optim(par=lambdamin,fn=PhiFL,method="CG")
  lambdamin=ROU$par
 }
  #plot( PhiL1[1:(lmax-1)],type="l")
  #return(rbind(c(lambdamin,PhiMin),c(ROU$par,ROU$value)))
  return(c(lambdamin,PhiMin))
}



PredS1=rep(0,NTA)
PredS1=TrajA

PredS2=rep(0,NTA)
PredS2=TrajA

Lambda1=array(rep(0,n*(NTE+NTP)),dim=c(NTE+NTP,n))
Lambda1[ITP-min(ITE,ITP)+1,]=diag(1,n)
Lambda2=Lambda1


lmax<-300
nlr<-(10/(10+1:lmax))^0.7
for (j in 1:NTE)
{
  st.time<-Sys.time()
  vt0=ITE[j]
  if (j==1){
    lambda0<-runif(n)
  } else
  {
    lambda0<-Lambda1[ITE[j-1]-min(ITE,ITP)+1,]
  }
  REZ1<-LambdaMin(type1=1,ti=vt0,nl=nlr,lmax=lmax,gm1=0.9,lambda0)
  PredS1[vt0]=REZ1[1:n]%*%TrajA[ITP]
  Lambda1[vt0-min(ITE,ITP)+1,]=REZ1[1:n]
  #if (j==1){
    lambda0<-REZ1[1:n]
  #}  else
  #{
   # lambda0<-Lambda2[ITE[j-1]-min(ITE,ITP)+1,]
  #}
  REZ2<-LambdaMin(type1=4,ti=vt0,nl=nlr,lmax=lmax,gm1=5,lambda0)
  PredS2[vt0]=REZ2[1:n]%*%TrajA[ITP]
  Lambda2[vt0-min(ITE,ITP)+1,]=REZ2[1:n]
  #PredS2[vt0]=REZ[2,1:n]%*%TrajA[ITP]
  ##TrajA[vt0]
  #Lambda2[vt0-min(ITE,ITP)+1,]=REZ[2,1:n]
  print(j)
  Sys.time()-st.time
}

warnings()
plot(TA[NT0:NTA], TrajA[NT0:NTA], type = "l", col = "black")
points(TA[NT0:NTA], PredS1[NT0:NTA], type = "l", col = "red")
#points(TA[NT0:NTA], PredS2[NT0:NTA], type = "l", col = "darkgreen")
points(TP, TrajA[ITP],  col = "blue")


plot(TA[NT0:NTA], TrajA[NT0:NTA], type = "l", col = "black")
points(TA[NT0:NTA], PredS2[NT0:NTA], type = "l", col = "red")
#points(TA[NT0:NTA], PredS2[NT0:NTA], type = "l", col = "darkgreen")
points(TP, TrajA[ITP],  col = "blue")



#writing results
write.csv(cbind(sort(c(ITE,ITP)),Lambda1),file=paste("Lambda-U",ModelPred,filename2,sep = "", collapse = NULL))
write.csv(cbind(sort(c(ITE,ITP)),Lambda2),file=paste("Lambda-C",ModelPred,filename2,sep = "", collapse = NULL))

  }
}

plot(TA, TrajA, type = "l", col = "black")


ks.test(F1(Traj[]),runif(1500))


qqplot(y=F1(TrajA[]),x=seq(0,1,1/1500))
hist(F0(TrajA[]))
hist(TrajA,freq = FALSE)
plot(p0(sort(TrajA)))

qqplot(y=F0(Traj[]),x=seq(0,1,1/1500))
qqplot(y=F0(rnorm(1500)),x=seq(0,1,1/1500))
qqplot(y=F1(Traj[20*(1:75)]),x=seq(0,1,1/75))
plot(Traj[20*(1:75)])
abline(0,1)
plot(F0(rno))

library(MASS)
fitdistr(Traj, "Levy")
#  F1<-function(x) return(pnorm(x,mean=-0.48,sd=1.11))
F1sigma<-function(s)
  {
  F1xhat<-2-2*pnorm(sqrt(s/TrajA))
  return(mean(1/3+F1xhat^2-F1xhat)-total.dist(F1xhat)/(length(TrajA))^2)
}

F1<-function(x) return(2-2*pnorm(sqrt(0.836/x)))
0.8359422

nlm(f=F1sigma,p=0.0835)

F0<-function(x) return()



##Plots
PredS1Plot=read.csv("InterPresS1n10h002.csv")
PredS2Plot=read.csv("InnerPresS2n10h002.csv")




plot(TA, TrajA, type = "l", col = "black")
points(T0, Traj, type = "l", col = "blue")
points(TA[NT0:NTA], PredS1[NT0:NTA], type = "l", col = "red")
points(TP, TrajA[ITP], type = "o", col = "blue")


plot(TA, TrajA, type = "l", col = "black")
points(T0, Traj, type = "l", col = "blue")
points(TA[NT0:NTA], PredS2[NT0:NTA], type = "l", col = "red")
points(TP, TrajA[ITP], type = "o", col = "blue")








plot(TA[NT0:NTA], TrajA[NT0:NTA], type = "l", col = "black")
points(TA[NT0:NTA], PredS2[NT0:NTA], type = "l", col = "red")
points(TP, TrajA[ITP],  col = "blue")




Krigs=read.csv("InterRezulth002.csv",header = TRUE)




#Krigs=read.csv("ExtraRezulth002.csv",header = TRUE)

plot(TA[NT0:NTA], PredS2Plot[NT0:NTA,2], type = "l", col = "black")
points(TA[NT0:NTA], PredS2Plot[NT0:NTA,3], type = "l", col = "red")
points(TA[(NT0+1):NTA], Krigs$Xok, type = "l", col = "darkgreen")
points(TP,  PredS2Plot[ITP,2],  col = "blue")

plot(TA[NT0:NTA], PredS1Plot[NT0:NTA,2], type = "l", col = "black")
points(TA[NT0:NTA], PredS1Plot[NT0:NTA,3], type = "l", col = "red")
points(TA[(NT0+1):NTA], Krigs$Xkm, type = "l", col = "darkgreen")
points(TP,  PredS1Plot[ITP,2],  col = "blue")



dist((1:3),upper=FALSE)

cumsum(cummax(c(3:1, 2:0, 4:2)))