function [B, a] = direct(xx, y, h, nd)

%Searching the effective dimension reduction subspace of model
%  			y = g(B^TX) + e
%Useage: directions(x, y, h, d, g)
%Input: 
%     x --- expanaltory variables
%     y --- response
%     h --- bandwidth
%     d --- working dimension of the space
%Output:
%     Directions B
%Reference: Y. Xia, H. Tong, W.K. Li and L.X. Zhu, 
% "adaptive estimation of the effictive dimension space",(2001)  
% ------------------------------------------------------------
% Example
%	n = 200;
%	x = randn(n,4);
%	beta1 = [1 2 3 0]';
%	beta1 = beta1/sqrt(beta1'*beta1);
%	beta2 = [-2 1 0 1]';
%	beta2 = beta2/sqrt(beta2'*beta2);
%	y = (x*beta1).^2 + x*beta2 + 0.2*randn(n,1);
%	B = OPG(x, y, 0.5, 2)
%  % estimation errors
%	B0 = [beta1 beta2];
%	error = B'*(eye(4)-B0*B0')*B
[n,p] = size(xx);
[ny, py] = size(y); 

[n,p] = size(xx);
[ny, py] = size(y); 
mm = mean(xx);
xx = xx - repmat(mm,n,1);
ss = inv(xx'*xx/n)^0.5;
xx = xx*ss;


if (ny ~= n)
   disp('Error: matrix of x and y dont match')
   return
end
if (p <= 1)
   disp('Error: please check your matrix x')
   return
end
if (py > 1)
   disp('Error: please check your matrix y')
   return
end
if (sum(abs(mean(xx,1))) > 0) | (abs(mean(std(xx,1))-1)>0) 
%   disp('Warning: covariates must be standardized')
end;   
onep = ones(p,1);
onen = ones(n,1);
B = eye(p);

for iter = 1:5;
	ab = ones(p,n);
    a = y;
	for i = 1:n;
   	xi = xx - repmat(xx(i,:),n,1);
   	kernel = exp(-sum((xi*B).^2,2)/(2*h*h));
   	onexi = [onen xi];
   	xk = onexi.*repmat(kernel, 1, p+1);
   	abi = inv(xk'*onexi+eye(p+1)*0.0001)*xk'*y;
   	ab(:,i) = abi(2:p+1);
   	a(i) =  abi(1);
	end;
	[B0 D] = eig(ab*ab');
	[D I] = sort(diag(D));
	B = B0;
	for k = 1:p;
   	B0(:,k) = B(:,I(k));
	end
	for k = 1:nd;
   	  B(:,k) = B0(:,p-k+1);
	end;  
   B = [B(:,1:nd)];
end;


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