# This code finds the process H_K which is the invelope of Y_K the (k-1)-fold integral of
#a two sided Brownian Motion + t^{2K} when K is even (K >=2).
# m is the precision of the Brownian Motion Approximation by Schauder functions.

#For K=6, c=4, p1=20,p2=16366
#For K=6, c=6, p1=30,p2=24548 (we can't go beyond ~ 2*10^{-7})


IterativeSHk <- function(K=6,C=4,m=11,eps=10^{-7},p=20,p1=10,p2=16365){

grid <- seq(-C,C,2^{-m})
#IntBr <- Intbrownk8c6m11[4097:20481,]
IntBr <- intbrownk6c4m11
IntBr <- t(IntBr)
Mat0 <- matrix(0,nrow=2*K + p, ncol=2*K+p)
L.g <- length(grid)
# 1 is the location of the successive derivative of Y at -C, L.g is that of ...of at C.
Yd <- rbind(IntBr[,1],IntBr[,L.g])

# Select only the even derivatives of Y at -C and C.

Yd <- Yd[,seq(2,K,2)]

# this vector stores in the first row: Y^{(k-1)}(-c),Y^{(k-2)}(-c),...,Y(-c)
# and Y^{(k-1)}(c),Y^{(k-2)}(c),...,Y(c).

S0 <- c(-C,C)	
Alpha0 <- StartingSplineHk(K,C,Yd)
Coef0 <- Alpha0[1]
H <- EvaluateGrid(K,Alpha0,S0,grid)
Diff <- H - IntBr[K,]

# For later, we need to have the initial conditions in the "right form" (as it is required for ComputeSplineHk)
# hence, we need to reverse the components of Yd so that we start by Y(-/+ c) and finish by Y^{(k-2)}(-/+ c).

Yd.rev <- Yd
Yd.rev[1,] <- rev(Yd[1,])
Yd.rev[2,] <- rev(Yd[2,])

# Check whether H >= Y.

min.Diff <- min(Diff[p1:p2])
print(min.Diff)
Count <- 0
while(min.Diff < -eps){	
Count <- Count + 1
cat("Main Loup numb = ", Count, "\n")	
Diff.sort <- rank(Diff[p1:p2])
min.rank <- min(Diff.sort)
min.pos <- match(min.rank,Diff.sort)
thetamin <- grid[p1:p2][min.pos] # locate t*.
valmin <- Diff[p1:p2][min.pos]

print(c(thetamin,valmin))
#rm(Diff)
#rm(H)

# Compute the new spline for the new set of knots.

S <- c(S0,thetamin)
S <- sort(S)
print(S)
#locate the knots in the grid.

positions <- match(S,grid)	


# Here, we need to create the vector Y of boundary conditions as it should be given in ComputeSpline2.

Y <- InitialCondHk(K,C,Yd.rev,IntBr,positions)
p <- length(S)-2
Alpha <- ComputeSplineHk(K=K,Y=Y,S=S,Mat0)
Alpha <- as.numeric(Alpha)
Coef <- c(Alpha[1],Alpha[(2*K+1):(2*K+p)])
Coef <- cumsum(Coef)
min.C <- min(diff(Coef))
count <- 0
	
	   while(min.C < 0){
		count <- count+1
		cat("Sub loup numb = ",count," of the main loop numb=", Count, "\n")
	   	index <- IndexFuncHk(S0=S0,S=S,Coef0=Coef0,Coef=Coef)
		S <- S[-index]
		p <- length(S)-2
		positions <- match(S,grid)	
		Y <- InitialCondHk(K,C,Yd.rev,IntBr,positions)
		Alpha <- ComputeSplineHk(K=K,Y=Y,S=S,Mat0)
		Alpha <- as.numeric(Alpha)
       Coef <- c(Alpha[1],Alpha[(2*K+1):(2*K+p)])
		Coef <- cumsum(Coef)
		min.C <- min(diff(Coef))
	   }#while min.C < 0 
	  
H <- EvaluateGrid(K,Alpha,S,grid)
Diff <- H - IntBr[K,]
min.Diff <- min(Diff[p1:p2])

S0 <- S
Coef0 <- Coef

}#while min.Diff < -eps
print(Alpha)
print(positions)
#Mat.H <- H
#for(d in 1:(2*K-2)){
#Mat.H <- rbind(Mat.H,EvaluateGridDer(K,Alpha,S,grid,d))
#}
#Mat.H
VecK <- EvaluateGridDerOdd(K,Alpha,S,grid,K)
VecK[8193]

}#end of the function