function [gamma SEgamma B SEB G seG L seL] = dmave(x, y, nd);

% estimation of multiple-index model
%  Y = \gamma^\top X + G(B^\top X) + \varepsilon
%
% Input:
%     X: nxp, covariates
%     Y: nx1, response
%     nd: dimension of nonlinear part
%
% Output:
%    gamma: linear coefficient, SEgamma: SE of gamma
%    B    : dimension reduction direction in the nonlinear part.
%    SEB  : SE of B
%    G    : the link function
%    seG  : SE of G



[n,p] = size(x);
onen = ones(n,1);

mm = mean(x);
x = x - repmat(mm,n,1);
x0 = x;
ss = inv(x'*x/n)^0.5;
x = x*ss;

m = p;
B = eye(p);
Ba = B;
BI = Ba;
B = B(:,1:m);
noip0 = 0;
noip1 = 1;
iterstop = 0;
Btmp = B;
rige = std(y)*mean(std(x,[],1));
%y = (y-mean(y))/std(y);
%yfit = y;

[beta B] = Mimd2(x,y,nd);
m = nd;

x0 = x;
xc = x;
enn = 1/n/n/n;

niter = floor(p*3/2)*2;
for iter = 1:niter;
    x = xc;
    if iter >= p;
        x = x0;
    end
  	adj = p^(1/iter)*n^(1/(nd+4)-1/(m+4));
    h = p^(1/(iter+1))*mean(std(x*Ba,[],1))/n^(1/(m+4));
%    h = cvm(x*B, y-x*beta);

    h2 = 2*h*h*adj;
    
    
%    ABI = zeros(m, n*n);
    onexi = ones(n, m+1);
    for iiter = 1:4;
	dd = 0;
   	dc = 0;
        for j = 1:n;
        xij = x - repmat(x(j,:),n,1);
        sxij = sum((xij*Ba).^2,2) + sum(xij.^2,2)/1.5^iter;
        sd = sort(sxij);
      	ker = exp(-sxij/h2);  
      	ker = ker/sum(ker); 
        f = 1;
        %        f = sum(exp(-sum((xij*Ba).^2,2)/h2)); 

        rker = repmat(ker, 1, p+1);
   		onexi(:,1:m) = xij*B;
   		xk = (onexi.*rker(:, 1:m+1))';
   		abi = inv(xk*onexi+eye(m+1)/n)*(xk*(y-x*beta));
      
        xdij = [x  kron(abi(1:m)', xij)];
        kxdij = (repmat(ker, 1, (m+1)*p).*xdij)';
        dd = dd + kxdij*xdij*f;
        dc = dc + kxdij*(y-abi(m+1))*f;
%        yfit(j) = abi(m+1);
        
        end

       BETA = inv(dd + eye(length(dc))*enn)*dc;
       
   	   B0 = reshape(BETA(p+1:p+p*m), p, m); %*ABI;
       [B, R] = eig(B0*B0');
       B = B(:,p-m+1:p);
       Ba = B;
       if (max(svd(B*B'-BI*BI')) < 0.001);
           break
       end
       BI = B;
       beta = (eye(p)-B*B')*BETA(1:p);       
    end

   if (max(svd(B*B'-BI*BI')) < 0.001)*(iter>p+3);
       'ok'
       break
   end
   BI = Ba;   
   
end

yfit = y;
xB = x*B;
dd = 0;
g = y;
df = y;
h = h*1.5;
h2 = 2*h*h;
for j = 1:n;
      xBi = (xB - repmat(xB(j,:), n, 1))/h;      
      sxij = sum(xBi.^2,2);
%      sd = sort(sxij);
%      h22 = max(h2, sd(m+1));
      ker = exp(-sxij/h2)/(2*pi)^(nd/2)/h^nd;  
   	  onexi = [xBi onen];
   	  xk = (onexi.*repmat(ker, 1, m+1))';
      invxx = inv(xk*onexi+eye(m+1)*enn);
   	  abi = invxx*(xk*(y-x*beta));
      vi = invxx*(xk*x);
      fd(j,:) = abi(1:m)/h;
      vB(j,:) = vi(m+1,:)/h;
      g(j) = abi(m+1);
      df(j) = mean(ker);
      yfit(j) = abi(m+1);
%      fit0 = ker'*(y-x*beta)/(sum(ker)+1.0e-5);
%      ker(j) = 0;
%      fit1 = ker'*(y-x*beta)/(sum(ker)+1.0e-5);
%      yfit(j) = (fit0 + fit1)/2;
end
s2 = mean(df.*(y- x*beta -yfit).^2)/mean(df)

G = g;
seG = sqrt(0.2821*s2./(n*h^nd*df));



dd = x-vB;
dd0 = dd;
for i = 1:m
    dd = [dd repmat(fd(:,i),1,p).*dd0];
end
%BT = eye((m+1)*p);
%BT(1:p, 1:p) = eye(p) - B*B';
I = (df>2/(2*pi)^(nd/2)/h^nd/n);
I = df;
ddf = dd.*repmat(I, 1,size(dd,2));


dd2 = ddf'*ddf;
dd0 = dd'*ddf;


sss = kron(eye(m+1), ss);
[c d] = eig(dd0);
ddiag = diag(d);
ds = sort(ddiag);
I = find(ddiag<= ds(m*(m+1)));
ddiag(I) = 1;
ddiag = 1./ddiag;
ddiag(I) = 0;
invdd = c*diag(ddiag)*c';


S = sss*invdd*dd2*invdd*sss*s2;
SEa = sqrt(diag(S(1:p,1:p)));
alpha = ss*beta;


%%%%%%%%%%%%%%%%%%

df = y;
xa = x0*alpha;
L = y; df = y;
h = std(xa)/n^0.2;
h2 = 2*h*h;
for j = 1:n;
      xi = xa - xa(j);
      sxij = sum(xi.^2,2);
      ker = exp(-sxij/h2)/sqrt(2*pi)/h;  
      ker1 = ker.*xi;
      df(j) = mean(ker);
      c = (ker*(xi'*ker1) - ker1*sum(ker1));
      L(j) = c'*(y-g)/(sum(c)+1/n^2);
      df(j) = mean(ker);
end

seL = sqrt(0.2821*s2./(n*h*df));


%%%%%%%%%%%%%%%%%

B = ss*B;
SEb = B;
for i = 1:size(B,2);
    d = sqrt(B(:,i)'*B(:,i));
    B(:,i) = B(:,i)/d;
    A = S(i*p+1:(i+1)*p, i*p+1:(i+1)*p)/d^2;
    SEb(:,i) = sqrt(diag(A));
end


s = inv(B'*B)^0.5;
B1 = B*s;
A = (eye(p) - B1*B1');
alpha = A*alpha;
gamma = alpha;

C = [A  (-B1*kron(eye(m), alpha'))
    zeros(m*p, p) eye(m*p)];
S = C*S*C';

SEgamma = sqrt(diag(S(1:p,1:p)));