LV_weight_partial<-function(V,Xl,dtemp,h,q,Y,grad=T){
  N<-length(Y)
  if(is.vector(V)){k<-1}
  else{k<-length(V[1,])}
  Xlv<-Xl%*%V
  d<-dtemp-((Xlv^2)%*%rep(1,k))
  w<-exp(-0.5*(d/h)^2)
  w<-matrix(w,N,q)
  wn<-apply(w,2,sum)-rep(1,q)#new
  w<-apply(w,2,column_normalize)
  mY<-t(w)%*%Y
  sig<-t(w)%*%(Y^2)-(mY)^2
  W<-(kronecker(t(wn),rep(1,N)))##new
  if(grad==T){
    grad<-mat.or.vec(dim,k)
    tmp1<-(kronecker(sig,rep(1,N))-(as.vector(kronecker(rep(1,q),Y))-kronecker(mY,rep(1,N)))^2)
    if(k==1){
      grad_d<- -2*Xl*as.vector(Xlv)
      grad<-(1/h^2)*(1/sum(wn))*t(grad_d*as.vector(d)*as.vector(w)*as.vector(W))%*%tmp1 #new
      # wn_grad<-(-1/h^2)*t(grad_d*as.vector(d)*as.vector(w))%*%kronecker(diag(rep(1,q)),rep(1,N))
      # grad<- wn_grad%*%(sig-rep(var1[2],q))/(sum(wn))+grad
    }
    else{
      for (j in 1:(k)){
        grad_d<- -2*Xl*as.vector(Xlv[,j])
        grad[,j]<- (1/h^2)*(1/sum(wn))*t(grad_d*as.vector(d)*as.vector(w)*as.vector(W))%*%tmp1#new
        # wn_grad<-(-1/h^2)*t(grad_d*as.vector(d)*as.vector(w))%*%kronecker(diag(rep(1,q)),rep(1,N))
        # grad[,j]<- wn_grad%*%(sig-rep(var1[2],q))/(sum(wn))+grad
      }
    }
    ret<-list(t(wn)%*%sig/sum(wn),sig,grad)#new
    names(ret)<-c('var','sig','grad')
  }
  else{
    ret<-list(t(wn)%*%sig/sum(wn),sig)#new
    names(ret)<-c('var','sig')
  }
  
  return(ret)
}
################
stiefl_weight_partial_opt<-function(dat,h=NULL,k,k0=30,p=1,maxit=50,nObs=sqrt(length(dat[,1])),lambda_0=1,tol=10^(-3),sclack_para=0){
  Y<-dat[,1]
  X<-dat[,-1]
  N<-length(Y)
  dim<-length(X[1,])
  if(p<1){
    S<-est_varmat(X)
    tmp1<-q_ind(X,S,p)
    q<-tmp1$q
    ind<-tmp1$ind
  }
  else{
    q<-N
    ind<-1:N
  }
  Xl<-(kronecker(rep(1,q),X)-kronecker(X[ind,],rep(1,N)))
  dtemp<-apply(Xl,1,norm2)
  if(is.null(h)){
    S<-est_varmat(X)
    tr<-var_tr(S)
    h<-choose_h_2(dim,k,N,nObs,tr)
  }
  best<-exp(10000)
  Vend<-mat.or.vec(dim,k)
  sig<-mat.or.vec(q,1)
  for(u in 1:k0){
    Vnew<-Vold<-stiefl_startval(dim,k)
    #print(Vold)
    #print(LV(Vold,Xl,dtemp,h,q,Y)$var)
    Lnew<-Lold<-exp(10000)
    lambda<-lambda_0
    err<-10
    count<-0
    count2<-0
    while(err>tol&count<maxit){
      #print(Vold)
      tmp2<-LV_weight_partial(Vold,Xl,dtemp,h,q,Y)
      G<-tmp2$grad
      Lold<-tmp2$var
      W<-G%*%t(Vold)-Vold%*%t(G)
      stepsize<-lambda#/(2*sqrt(count+1))
      Vnew<-solve(diag(1,dim)+stepsize*W)%*%(diag(1,dim)-stepsize*W)%*%Vold
      # print(Vnew)
      tmp3<-LV_weight_partial(Vnew,Xl,dtemp,h,q,Y,grad=F)
      Lnew<-tmp3$var
      err<-sqrt(sum((Vold%*%t(Vold)-Vnew%*%t(Vnew))^2))/sqrt(2*k)#sqrt(sum(tmp3$grad^2))/(dim*k)#
      #print(err)
      if(((Lnew-Lold)/Lold) > sclack_para){#/(count+1)^(0.5)
        lambda=lambda/2
        err<-10
        count2<-count2+1
        count<-count-1
        Vnew<-Vold #!!!!!
        
      }
      Vold<-Vnew
      count<-count+1
      #print(count)
    }
    if(best>Lnew){
      best<-Lnew
      Vend<-Vnew
      sig<-tmp3$sig
    }
  }
  ret<-list(Vend,best,sig,count,h,count2)
  names(ret)<-c('est_base','var','aov_dat','count','h','count2') 
  return(ret)
}

#################MAVE, OPG, rMAVE, rOPG from Bing Li book
opg=function(x,y,d){
  p=dim(x)[2];
  n=dim(x)[1] 
  c0=2.34;
  p0=max(p,3);
  rn=n^(-1/(2*(p0+6)));
  h=c0*n^(-(1/(p0+6)))
  sig=diag(var(x));
  x=apply(x,2,standvec)
  kmat=kern(x,h);
  bmat=numeric() 
  for(i in 1:dim(x)[1]){ 
    wi=kmat[,i];
    xi=cbind(1,t(t(x)-x[i,]))
    bmat=cbind(bmat,wls(xi,y,wi)$b)}
  beta=eigen(bmat%*%t(bmat))$vectors[,1:d] 
  return(diag(sig^(-1/2))%*%beta)
}
#######################
stiefl_opt_momentum<-function(dat,h=NULL,k,k0=30,p=1,maxit=50,nObs=sqrt(length(dat[,1])),lambda_0=1,tol=10^(-3),sclack_para=0,momentum_para=0.8){
  Y<-dat[,1]
  X<-dat[,-1]
  N<-length(Y)
  dim<-length(X[1,])
  if(p<1){
    S<-est_varmat(X)
    tmp1<-q_ind(X,S,p)
    q<-tmp1$q
    ind<-tmp1$ind
  }
  else{
    q<-N
    ind<-1:N
  }
  Xl<-(kronecker(rep(1,q),X)-kronecker(X[ind,],rep(1,N)))
  dtemp<-apply(Xl,1,norm2)
  if(is.null(h)){
    S<-est_varmat(X)
    tr<-var_tr(S)
    h<-choose_h_2(dim,k,N,nObs,tr)
  }
  best<-exp(10000)
  Vend<-mat.or.vec(dim,k)
  sig<-mat.or.vec(q,1)
  for(u in 1:k0){
    Vold<-stiefl_startval(dim,k)
    #print(Vold)
    #print(LV(Vold,Xl,dtemp,h,q,Y)$var)
    Lnew<-Lold<-exp(10000)
    lambda<-lambda_0
    err<-10
    count<-0
    count2<-0
    Lnew<-LV(Vold,Xl,dtemp,h,q,Y)$var
    #print(Lnew)
    if(best>Lnew){
      best<-Lnew
      Vend<-Vold
      #sig<-tmp3$sig
    }
  }
  Vnew<-Vold<-Vend
  
  G<-matrix(rep(0,dim*k),dim,k)
  while(err>tol&count<maxit){
    #print(Vold)
    tmp2<-LV(Vold,Xl,dtemp,h,q,Y)
    #G<-tmp2$grad
    G<-(1-momentum_para)*G +  momentum_para*tmp2$grad
    Lold<-tmp2$var
    W<-G%*%t(Vold)-Vold%*%t(G)
    stepsize<-lambda#/(2*sqrt(count+1))
    Vnew<-solve(diag(1,dim)+stepsize*W)%*%(diag(1,dim)-stepsize*W)%*%Vold
    # print(Vnew)
    tmp3<-LV(Vnew,Xl,dtemp,h,q,Y,grad=F)
    Lnew<-tmp3$var
    err<-sqrt(sum((Vold%*%t(Vold)-Vnew%*%t(Vnew))^2))/sqrt(2*k)#sqrt(sum(tmp3$grad^2))/(dim*k)#
    #print(err)
    if(((Lnew-Lold)/Lold) > sclack_para){#/(count+1)^(0.5)
      lambda=lambda/2
      err<-10
      count2<-count2+1
      count<-count-1
      Vnew<-Vold #!!!!!
      Lnew<-Lold
      
    }
    Vold<-Vnew
    count<-count+1
    
    #print(count)
  }
  
  ret<-list(Vnew,Lnew,count,h,count2)
  names(ret)<-c('est_base','var','count','h','count2') 
  return(ret)
}
######################
wls=function(x,y,w){ 
  n=dim(x)[1];
  p=dim(x)[2]-1 
  out=c(solve(t(x*w)%*%x/n)%*%apply(x*y*w,2,mean))
  return(list(a=out[1],b=out[2:(p+1)]))
}
#################
kern=function(x,h){
  x=as.matrix(x);
  n=dim(x)[1]
  k2=x%*%t(x);
  k1=t(matrix(diag(k2),n,n));
  k3=t(k1);
  k=k1-2*k2+k3 
  return(exp(-(1/(2*h^2))*(k1-2*k2+k3)))
}
###############
standvec=function(x) return((x-mean(x))/sd(x))
##############
mave2=function(x,y,h,d,nit){ 
  sig=diag(var(x));
  n=dim(x)[1];
  p=dim(x)[2] 
  x=apply(x,2,standvec);
  beta=opg(x,y,d);#beta=opg(x,y,h,d);
  kermat=kern(x,h) 
  for(iit in 1:nit){ 
    b=numeric();
    a=numeric();
    for(i in 1:n){ 
      wi=kermat[,i]/(apply(kermat,2,mean)[i])
      ui=cbind(1,t(t(x)-x[i,])%*%beta) 
      out=wls(ui,y,wi);
      a=c(a,out$a);b=cbind(b,out$b)}
    out=0;out1=0;
    for(i in 1:n){ 
      xi=kronecker(t(t(x)-x[i,]),t(b[,i]))
      yi=y-a[i];
      wi=kermat[,i]/apply(kermat,2,mean)[i]
      out=out+apply(xi*yi*wi,2,mean) 
      out1=out1+t(xi*wi)%*%xi/n} 
    beta=t(matrix(solve(out1)%*%out,d,p))
  } 
  return(diag(sig^(-1/2))%*%beta)
}
######################
rmave=function(x,y,d,nit){ 
  sig=diag(var(x));
  n=dim(x)[1];
  p=dim(x)[2]
  x=apply(x,2,standvec) 
  c0=2.34;
  p0=max(p,3);
  h=c0*n^(-(1/(p0+6)));
  rn=n^(-1/(2*(p0+6))) 
  beta=opg(x,y,d) 
  for(iit in 1:nit){ 
    kermat=kern(x%*%beta,h);
    mkermat=apply(kermat,2,mean) 
    b=numeric();a=numeric() 
    for(i in 1:n){ 
      wi=kermat[,i]/mkermat[i];
      ui=cbind(1,t(t(x)-x[i,])%*%beta) 
      out=wls(ui,y,wi);
      a=c(a,out$a);b=cbind(b,out$b)
    } 
    out=0;
    out1=0 
    for(i in 1:n) { 
      xi=kronecker(t(t(x)-x[i,]),t(b[,i]));
      yi=y-a[i] 
      wi=kermat[,i]/mkermat[i] 
      out=out+apply(xi*yi*wi,2,mean) 
      out1=out1+t(xi*wi)%*%xi/n} 
    beta=t(matrix(solve(out1)%*%out,d,p)) 
    h=max(rn*h,c0*n^((-1/(d+4)))) 
  } 
  return(diag(sig^(-1/2))%*%beta)
}
#########################
ropg=function(x,y,d,nit){ 
  sig=diag(var(x));
  x=apply(x,2,standvec);
  p=dim(x)[2];
  n=dim(x)[1]
  c0=2.34;
  p0=max(p,3);
  rn=n^(-1/(2*(p0+6)));
  h=c0*n^(-(1/(p0+6))) 
  beta=diag(p) 
  for(iit in 1:nit){ 
    kmat=kern(x%*%beta,h);
    bmat=numeric() 
    for(i in 1:dim(x)[1]){ 
      wi=kmat[,i];
      xi=cbind(1,t(t(x)-x[i,])) 
      bmat=cbind(bmat,wls(xi,y,wi)$b)
    } 
    beta=eigen(bmat%*%t(bmat))$vectors[,1:d] 
    h=max(rn*h,c0*n^((-1/(d+4)))) 
  }
  
  beta.final=diag(sig^(-1/2))%*%beta 
  return(beta.final) 
}