%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This code is to compute SCMI predictive region using the kNN local
linear smoothing
approach
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matrixprod2<-function (A) {

result<-array(0, dim=c(nrow(A),nrow(A), ncol(A)))

 for( k in 1:ncol(A))
   {
                       result[,,k]<- outer(A[,k], A[,k], "*"
)}
                       return(result)
                       }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
funopare.Shortset.lcv <- function(Response,
                                  CURVES,
                                  ...,
                                  alpha=0.9,
                                  Knearest=NULL,
                                  PRED=CURVES,
                                  kind.of.kernel = "quadratic",
                                  semimetric = "deriv")
{
        sm <- get(paste("semimetric.", semimetric, sep=""))
        kernel <- get(kind.of.kernel)
           n <- nrow(CURVES)
           SMLEARN1 <- sm(CURVES,...)

 zero<-matrix(0,n,length(Theta1))


        SMLEARN2NEWtheta<-matrix(0,n,n)
for (i in 1:n)
                 SMLEARN2NEWtheta[,i] <-
                 ((sm(CURVES-CURVES[i,],zero[i,],...))^2+(sm(zero[1:n,],Theta1,...))^2-(sm(CURVES-CURVES[i,],Theta1,...))^2)/2

Wij<- array(0,dim=c(n,n,n))

for (i in 1:n)
Wij[,,i]<-SMLEARN2NEWtheta[,i]*outer(SMLEARN2NEWtheta[,i],SMLEARN2NEWtheta[,i],"-")




           SML12.SOR <- apply(SMLEARN1, 2, sort)
            Resp.range <- range(Response)
        Response.grid <- seq(from=Resp.range[1]*0.9,to=Resp.range[2]*1.1,
                             length=100)
        # RESPMETRIC[i,j]=yi-yj with i in Response.grid and j in LEARN1
        RESPMETRIC <- outer(Response.grid, Response, "-")

         lgrid <- length(Response.grid)

if(is.null(Knearest)){
        Knearest.min <- floor(n *0.05)
                Knearest.max <- floor(n * 0.25)
                step <- ceiling((Knearest.max - Knearest.min) / 10)
                Knearest <- seq(Knearest.min, Knearest.max, by = step)
       }
            lknearest <- length(Knearest)

BANDL12.CUR <- 0.5 * (SML12.SOR[Knearest,] +
SML12.SOR[Knearest+1,])

CV <- matrix(0, nrow = lknearest, ncol = n)


        QUANT <- matrix(0, nrow = lknearest, ncol = n)

 count1 <- 0
        for(kk in Knearest){
            count1 <- count1 + 1
            ARG <- t(t(SMLEARN1)/BANDL12.CUR[count1,])
            KERNEL.CURVES <- kernel(ARG)
        KERNEL.CURVES[KERNEL.CURVES<=0] <- 0
        KERNEL.CURVES[KERNEL.CURVES>=1] <- 0
        Denom <- apply(KERNEL.CURVES, 2, sum)

            IKERNEL.RESP  <- apply((RESPMETRIC), 1, sign)
 IKERNEL.RESP[ IKERNEL.RESP<=0] <- 0
            CDF.ESTIMATE <- (t(KERNEL.CURVES)/Denom) %*%
                                                  IKERNEL.RESP
            CDF.ALPHA <- CDF.ESTIMATE - 0.5
            Ind.quant <- apply(CDF.ALPHA < 0, 1,sum) + 1
            QUANT[count1,] <- Response.grid[Ind.quant]
            CV[count1,] <- abs(Response-QUANT[count1,])
        }


    Ind.knearest.opt <- apply(CV, 2, order)[1,]
IND.OPT <- cbind(Ind.knearest.opt,1:n)
    Bandwidth.opt<-BANDL12.CUR[IND.OPT]

 if(!missing(PRED)){
                SMLEARN2NEW <- sm(CURVES,...,PRED)
                Order.new <- apply(SMLEARN2NEW, 2, order)[1,]
                Bandwidth2 <- 0
                Ext.Short.predi.inter.g<-0
                  Ext.Short.predi.inter.d<-0


                n2 <- ncol(SMLEARN2NEW)

                                qunt.pridic<-matrix(0,n2,lgrid)
                                Leb<-matrix(0,lgrid, lgrid)
                for(k in 1:n2){
                 Bandwidth2<-Bandwidth.opt[Order.new[k]]
    KERNEL <-kernel((SMLEARN2NEW[,k])/Bandwidth2)
                KERNEL[KERNEL<0] <- 0
                KERNEL[KERNEL>=1] <- 0
                Denom1 <- sum(KERNEL)
                while(Denom1==0) {
                        Bandwidth2 <- 1.1 * Bandwidth2
                        KERNEL <- kernel((SMLEARN2NEW[,k])/Bandwidth2)
                        KERNEL[KERNEL<0] <- 0
                KERNEL[KERNEL>=1] <- 0
                Denom1 <- sum(KERNEL)
                }

               qunt.pridic[k,] <- (t(KERNEL)/Denom1)%*% (IKERNEL.RESP)

Dif.Fun.rep.pridic<- abs(outer(qunt.pridic[k,],qunt.pridic[k,],
"-"))

Dif.rep.pridic<-abs(outer(Response.grid,Response.grid, "-"))

Leb[Dif.Fun.rep.pridic>alpha]<-Dif.rep.pridic[Dif.Fun.rep.pridic>alpha]

Leb[Dif.Fun.rep.pridic<=alpha]<-max(Response.grid)-min(Response.grid)

ind.Ext.Inter<-order(t(Leb))[1]


Ext.Inter.g <-ceiling(ind.Ext.Inter/lgrid)
Ext.Inter.d<-ind.Ext.Inter-(ceiling(ind.Ext.Inter/lgrid)-1)*lgrid

Ext.Short.predi.inter.g[k]<-min(Response.grid[Ext.Inter.g],Response.grid[Ext.Inter.d])

Ext.Short.predi.inter.d[k]<-max(Response.grid[Ext.Inter.g],Response.grid[Ext.Inter.d])
}






                 Shoertset.Interval <-cbind(Ext.Short.predi.inter.g,Ext.Short.predi.inter.d)

                return(list(Shoert.Interval=Shoertset.Interval,
))
        }
        else{
                return(list(Response.values=Response
                            ))
        }

}