function theta = abc(x, y);

% y: response
% X = (x_1, ..., x_p)^T: covariates 
% model: y = g_1(b^T X) x_1 + ... + g_p(b^T X) x_p + e
%
% input----  x: nxp; y:nx1
% output---  theta: px1
%
% Example
%	 n = 100;
%	 p = 4;
%	 x =  2*(rand(n,p)-0.5);
%    
%    theta0 = [1 2 0 2]'/3;
%    x1 = x*theta0;
%    y = 3*exp(-x1.*x1) + 0.8*x1.*x(:,1) + 1.5*sin( pi*x1).*x(:,3)+ 0.5*randn(n,1);
%    
%    theta = sifc(x, y)
%
% Reference: 
%     Xia and Li, 1999, JASA, page 1275-1285
%     Fan, Yao and Cai, 2002, JRSS(B) page 57-80
%     Xia, Li and Tong, 2005, (to appear in Statistica Sinica)
%

xx = x;
[n, p] = size(xx);
zz = [ones(n,1) x];
[n, q] = size(zz);

h = mean(std(xx))/n^(1/(p+3));
theta = eye(p);

for ii = 1:10
h2 = 2*h*h;
for i = 1:n
    xi = xx - repmat(xx(i,:),n,1);
    k = exp( - sum((xi*theta).^2,2)/h2);
    xz = [];
    for j = 1:q
        xz = [xz xi.*repmat(zz(:,j),1,p)];
    end
    xz = [xz zz];
    zxk = xz.*repmat(k,1,(p+1)*q);
    ab(:,i) = inv(xz'*zxk+eye((p+1)*q)/sqrt(n))*(zxk'*y);
end
ab0 = [];
for j = 1:q;
    ab0 = [ab0 ab((j-1)*p+1:j*p,:)];
end
[D, v] = eig(ab0*ab0');
v = diag(v);
I = max(find(v == max(v)));
theta = D(:,I);
h = mean(std(xx*theta))/n^(1/(p+3));
end

f = zz;
g = zz;
niter = 1;
nref =  5;
b = theta;

for ii = 1:niter
	xx0 = xx*b/sqrt(b'*b);
	h0 = std(xx0)*n^(-.2);
   h0 = sifc_cv(0.5*h0, 4*h0,  xx0, zz, y);
	h2 = 2*h0*h0;
   b0 = b;
	Jp = ones(size(xx0,2),1);
	for iter = 1 : niter;
		mc0 = 0.0;
		md0 = 0.0;
        cve = 0.0;
		xx1 = xx*b;
		for j = 1 : n;
			zxx = xx0 - repmat(xx0(j,:),n,1);
			k = exp(- zxx.*zxx*Jp/h2);
			k = k/sum(k);
			zxt = [zz zz.*repmat( xx1 - xx1(j,:), 1, q)];
			kzxt = repmat(k,1,2*q).*zxt;
			ma = kzxt'*zxt + eye(2*q)*0.00001;
			mb = kzxt'*y;
            ab = inv(ma)*mb;
            f(j,:) = ab(1:q)';
            g(j,:) = ab(q+1:2*q)';
           
			zxx = xx - repmat(xx(j,:),n,1);
			zbx = repmat(zz*ab((q+1):(2*q)), 1, p).*zxx;
			kzbx = repmat(k, 1, p).*zbx;
			mc0 = mc0 + kzbx'*zbx;
			md0 = md0 + kzbx'*(y-zz*ab(1:q));
		end;
		b = inv(mc0+eye(p)/n)*md0;
		b = ( (b(2)>0 ) - ( b(2) < 0) )* b /sqrt(b'*b);
   end;	

end
   
b = b/sqrt(b'*b);
theta = b;

%************************************

function h0 = sifc_cv(hL, hU, xx, zz, y);

	[n, p] = size(xx);
    [n, q] = size(zz);
    f = zz;
   
    ss0 = 1.0e50;
	hD = (hU - hL);
	for ih = 1:10;
	ss = 0.0;
	h = hL + hD*ih/10;
	hm = 2*h*h;
    for i = 1 : n;
        	xy = xx - repmat(xx(i, :), n, 1);
        	k = exp(- sum(xy.*xy, 2)/hm );
         k(i) = 0.0;
      	k = repmat(k/( sum(k) +1.0e-10 ), 1, q);
         kzz = (k.*zz)';
         va = kzz*zz + 1.0E-4*eye(q);
         vb = kzz*y;
         a = inv(va)*vb;
	      ss = ss + (y(i) - zz(i,:)*a)^2;
	end;
	if ss < ss0;
      ss0 = ss;
      h0 = 1.5*h;  % 1.5 is the adjustment coefficient between NW and LL kernel smoothers
	end;
end;