function [B, cv] = dmave(x, y, m0);

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

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));

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

y = (y-mean(y))/std(y);
x = x./repmat(std(x,[],1),n,1);

yc = y;
x0 = x;
xc = x;

if 1 == 1
[a, yt] = sort(y);
[a, yt] = sort(yt);
ycc = yt/max(yt);
for i = 1:size(x,2);
    a = x(:,i);   
    [b, a] = sort(a);
    [b, a] = sort(a);
    xc(:,i) = a/std(a);
end
end


adj0 = 1;
niter = floor(p*3/2);
ABI = zeros(m, n*n); 
%ch = (sqrt(p)/n^(1/(p+4))*n^(2/(m0+4)))^(1/niter);
for iter = 1:niter;
    x = xc;
    if iter >= p;
        x = x0;
    end
  	adj = p^(1/iter)*n^(1/(m0+4)-1/(m+4));
  	hy = std(yc)/n^(1/(m0+4))*adj*adj0;
	ky = repmat(yc, 1, n);
    ky = [ycc/sqrt(2*pi)/hy*n^(1/iter) exp(-(ky-ky').^2/(2*hy*hy))/sqrt(2*pi)/hy];
    n1 = size(ky,2);
    h = p^(1/(iter+1))*mean(std(x*Ba,[],1))/n^(1/(m+4));
    h2 = 2*h*h*adj*adj0;
    
   
   	BB0 = 0;
   	for j = 1:n;
        xij = x - repmat(x(j,:),n,1);
        dxij = sum((xij*Ba).^2,2);
        h2i = h2;
        if length(find(dxij<h2))< max(5, sqrt(n)/2);
            h2i = h2*1.2;
        end
        ker0 = dxij + sum(xij.^2,2)/2^iter;
      	ker = exp(-ker0/h2i);      
   		onexi = [xij onen];
   		xk = (onexi.*repmat(ker, 1, p+1))';
   		abi = inv(xk*onexi+eye(p+1)/n)*xk*ky;
        BB0 = BB0 + abi(1:m,:)*abi(1:m,:)';
   	end
	   
       [B, R] = eig(BB0/n);
       R = sqrt(R);
       r1 = max(diag(R)');
       Ba = B*R*B'/r1;
       
end

[R, I] = sort(diag(R));

B = B(:,I);
B = B(:, p-m0+1:p);

B = ss*B;
for i = 1:size(B,2);
    B(:,i) = B(:,i)/sqrt(B(:,i)'*B(:,i));
end