# This function is equivalent to the Iterative (2k-1)th spline. The goal is, given 
# a k-monotone function f (the current iterate), to find the next iterate as a linear
#combination of a (2k-1)-th spline and the current iterate. Each time a support point
# is deleted from the set of knots. Few iterations are necessary till we get a spline
# P such that P^{(k)} is a k-monotone function.


LSESupReducAlgo <- function(K=3,X=X1000.3,prec=0.01,eps= 10^{-8},p1=1,p2=1){
	
	#theta0 <- (2*K-1)*max(X)
	grid <- round(seq(min(X)*p1,p2*K*max(X),prec),digits=6)
	M.alpha <- matrix(0,nrow=K-1,ncol=K-1) 
	M0 <- matrix(0,nrow=30,30)
	B0 <- rep(0,30)
	#grid <- round(seq(min(X),2*K*max(X),prec),digits=6)
	Rank <- rank(c(max(X),grid))[1]
   	theta0 <- grid[Rank]
	#theta0 <- grid[length(grid)]
   	print(theta0)
	C0 <- ((2*K-1)/(K*theta0^{K-1}))*mean((theta0-X)^{K-1})
	#print(C0)	
	S0 <- theta0
	Resmin <- FindMinFunc(X=X, S=S0,C=C0,K=K,prec=prec,grid)
	valmin <- Resmin[2]
	print(valmin)
	Count <- 0
	while(valmin < -eps){  
	Count <- Count + 1	
	cat("Main loup numb = ",Count,"\n")
	thetamin <- Resmin[1]
	print(c(thetamin,valmin))
 	S <- c(S0,thetamin)
	S <- sort(S)
	B <- LSEInitialCond(S=S,K=K,X=X,B0)
	C <- LSEComputeSpline(S=S,K=K,B=B,M.alpha,M0)
	C <- ((-1)^K * S^K * factorial((2*K-1))/factorial(K))*C
	print(S)
	print(C)
	min.C <- min(C)
	count <- 0
	   while(min.C < 0){
		count <- count+1
		cat("Sub loup numb = ",count," of the main loop numb=", Count, "\n")
	   	index <- IndexFunc(S0=S0,S=S,C0=C0,C=C)
		S <- S[-index]
		if(length(S)==1)
		C <- ((2*K-1)/(K*S^{K-1}))*mean((S-X[X <= S])^{K-1})
		else{
		B <- LSEInitialCond(S=S,K=K,X=X,B0)
		C <- LSEComputeSpline(S=S,K=K,B=B,M.alpha,M0)
       C <- ((-1)^K * S^K * factorial((2*K-1))/factorial(K))*C
		}
		min.C <- min(C)
	  
    
}# while(min.C < 0)		
			S0 <- S
          	C0 <- C
			Resmin <- FindMinFunc(X=X, S=S0,C=C0,K=K,prec=prec,grid)
			valmin <- Resmin[2]

       
    }# while(valmin < -eps) 

       Output <- cbind(S0,C0)
       Output


}