function [Bs, Cs, BL, CL] = SCA(x, z, m)

% A code for semiparametric canonical analysis, Xia, Y. (2006)
% 
% input: x: nxp
%        z: nxq
%        m: specified number of components
%
% output: [Bs, Cs]: nonlinear component coefficients
%           (x*Bs(:,k), z*Cs(:,k)) is the kth NCA components  
%         [BL, CL]: linear component coefficients
%           (x*BL(:,k), z*CL(:,k)) is the kth LCA components  
%
%
%% Example
% n = 50;
% X = randn(n, 4);
% C = [.707  0  .707   0  0
%         0 .6     0  .8  0
%      .707  0 -.707   0  0
%         0 .8     0 -.6  0
%         0  0     0   0  1];
%
% B = [ .5   .5
%      -.5   .5
%       .5   .5
%      -.5   .5];
%
% Z = ([X*B(:,1) (X*B(:,2)).^2 zeros(n,3)] + randn(n,5)*0.2)*inv(C);
%
% [Bs, Cs, BL, CL] = SCA(X, Z, 2)
%
%% another example is at the bottom of the code


[n, p] = size(x);
[n,q] = size(z);
m = min([q p m]);

z = z - repmat(mean(z), n, 1);
x = x - repmat(mean(x), n, 1);
sz = (inv(z'*z/n))^.5;
z = z*sz;
sx = (inv(x'*x/n))^.5;
x = x*sx;
[C B] = LCA(x,z);
BL = B(:,1:m); CL = C(:,1:m);

Bs = [BL BL*0];
Cs = [CL CL*0];

C = Bs(:,1:m); 
B = Cs(:,1:m);
for j = 1:10
    for im = 1:m
        [C(:,im), B(:,im), yk(:,im)] = sCCA(x,z, C(:,[1:im-1 im+1:m]));
    end
    if m == 1
        break
    end    
end
Bs = B; Cs = C;
Cs = inv(z'*z)*(z'*yk);
for im = 1:m
        Cs(:,im) = Cs(:,im)/sqrt(Cs(:,im)'*Cs(:,im));
end

%figure

for im = 1:m
    [xtheta, I] = sort(x*Bs(:,im));
    yi = z(I,:)*Cs(:,im);
    xp = min(xtheta) + (0:100)'/100*(max(xtheta)-min(xtheta));
    h = std(xtheta)/n^0.2;
    yp = xp;
    f = yp;
    ye = yi;
    for j = 1:n
        xi = xtheta - xtheta(j);
        k = exp( -xi.^2/(2*h*h))/sqrt(2*pi)/h;
        k(j) = 0;
        k1 = xi.*k;
        s0 =sum(k); s1 = sum(k1); s2 = sum(k1'*xi);
        ye(j) = (s2*k-s1*k1)'*yi/(s0*s2-s1*s1+1.0e-10);
    end
    
    for j = 0:100
        xi = xtheta - xp(j+1);
        k = exp( -xi.^2/(2*h*h))/sqrt(2*pi)/h;
        k1 = xi.*k;
        s0 =sum(k); s1 = sum(k1); s2 = sum(k1'*xi);
        yp(j+1) = (s2*k-s1*k1)'*yi/(s0*s2-s1*s1+1.0e-10);
        f(j+1) = s0/n;
    end
    v = sqrt(0.2821./(n*h*f)*mean((yi-ye).^2));
    
    subplot(3,3,im)
    cc = plot(x*Bs(:,im), z*Cs(:,im), '.');
    hold on
    set(cc, 'markersize', 13)
    cc = title([num2str(im) 'th pair of SCA components']);  
    set(cc, 'fontsize', 12)
    
    cc = plot(xp, yp, 'r-');
    set(cc, 'linewidth', 2)
    cc = plot(xp, yp+2*v, 'r--');
    cc = plot(xp, yp-2*v, 'r--');
end

% output 

Cs = sz*Cs;
Bs = sx*Bs;
BL1 = sx*CL;
CL = sz*BL;
BL = BL1;

for im = 1:m
    Cs(:,im) = Cs(:,im)/sqrt(Cs(:,im)'*Cs(:,im));
    Bs(:,im) = Bs(:,im)/sqrt(Bs(:,im)'*Bs(:,im));
    CL(:,im) = CL(:,im)/sqrt(CL(:,im)'*CL(:,im));
    BL(:,im) = BL(:,im)/sqrt(BL(:,im)'*BL(:,im));
end

return

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

function [C, B, yk] = sCCA(xx, zz, Ba)

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

onep = ones(p,1);
onen = ones(n,1);

B = eye(p);
C = eye(q);



[Q, v] = eig(Ba*Ba');
v = diag(v);
Q = Q(:,find(v>1.0e-4));
Q = eye(q) - Q*Q';

    xB = xx*B;
    h0 = mean(std(xB))/n^(1/(size(B,2)+4));
    h = h0;
    zzx = zz;
    for i = 1:n
        xBi = (xB - repmat(xB(i,:),n,1))/h;
        k0 = exp(-sum(xBi.^2,2)/2);
        zzx(i,:) = k0'*zz/(sum(k0)+1/n);
    end
    abz = zz - zzx;
    [C0 D] = eig(abz'*abz);
	[D I] = sort(diag(D));
	C = C0;
	for k = 1:q;
    	C0(:,k) = C(:,I(k));
    end
    C = C0(:,1);    
    C = Q*C;
    C = C/sqrt(C'*C);

for iter = 1:10;    
    y = zz*C;
    h = h0/std(y);
    qC = size(C,2);
    ab = ones(p,n*qC);
    h2 = 2*h*h;
    yk = y;
	for i = 1:n;
        xi = xx - repmat(xx(i,:),n,1);
        kernel = exp(-sum((xi*B).^2,2)/h2  );
        onexi = [onen xi];
        xk = onexi.*repmat(kernel, 1, p+1);
        abi = inv(xk'*onexi+eye(p+1)/n)*xk'*y;
        yk(i) = abi(1,1);
        ab(:,(i-1)*qC+1:i*qC) = abi(2:p+1,:);
	end;
	[B0 D] = eig(ab*ab');
	[D I] = sort(diag(D));
	B = B0;
	for k = 1:p;
   	B0(:,k) = B(:,I(k));
	end
   	B = B0(:,p);
    
    zz1 = [zz ones(n,1)];
    C = inv(zz1'*zz1)*zz1'*yk;
    C = C(1:q);
    C = Q*C;
    C = C/sqrt(C'*C);

end;

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

function [Wx, Wy] = LCA(X,Y)

X = X'; Y = Y';
z = [X; Y];
C = cov(z.');
sx = size(X,1);
sy = size(Y,1);
Cxx = C(1:sx, 1:sx) + 10^(-8)*eye(sx);
Cxy = C(1:sx, sx+1:sx+sy);
Cyx = Cxy';
Cyy = C(sx+1:sx+sy, sx+1:sx+sy) + 10^(-8)*eye(sy);
invCyy = inv(Cyy);

[Wx,r] = eig(inv(Cxx)*Cxy*invCyy*Cyx); 
r = sqrt(real(r));     

V = fliplr(Wx);		
r = flipud(diag(r));	
[r,I]= sort((real(r)));	
r = flipud(r);		
for j = 1:length(I)
  Wx(:,j) = V(:,I(j)); 
end
Wx = fliplr(Wx);	

Wy = invCyy*Cyx*Wx;     
Wy = Wy./repmat(sqrt(sum(abs(Wy).^2)),sy,1); 

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Example 2: salesperson's data (John & Wichern, 2002, pp.539)  %
%% To run the example, just copy and paste to the command window %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

zx=[93.0  96.0  97.8  9  12  9  20
   88.8  91.8  96.8  7  10  10  15
   95.0  100.3  99.0  8  12  9  26
   101.3  103.8  106.8  13  14  12  29
   102.0  107.8  103.0  10  15  12  32
   95.8  97.5  99.3  10  14  11  21
   95.5  99.5  99.0  9  12  9  25
   110.8  122.0  115.3  18  20  15  51
   102.8  108.3  103.8  10  17  13  31
   106.8  120.5  102.0  14  18  11  39
   103.3  109.8  104.0  12  17  12  32
   99.5  111.8  100.3  10  18  8  31
   103.5  112.5  107.0  16  17  11  34
   99.5  105.5  102.3  8  10  11  34
   100.0  107.0  102.8  13  10  8  34
   81.5  93.5  95.0  7  9  5  16
   101.3  105.3  102.8  11  12  11  32
   103.3  110.8  103.5  11  14  11  35
   95.3  104.3  103.0  5  14  13  30
   99.5  105.3  106.3  17  17  11  27
   88.5  95.3  95.8  10  12  7  15
   99.3  115.0  104.3  5  11  11  42
   87.5  92.5  95.8  9  9  7  16
   105.3  114.0  105.3  12  15  12  37
   107.0  121.0  109.0  16  19  12  39
   93.3  102.0  97.8  10  15  7  23
   106.8  118.0  107.3  14  16  12  39
   106.8  120.0  104.8  10  16  11  49
   92.3  90.8  99.8  8  10  13  17
   106.3  121.0  104.5  9  17  11  44
   106.0  119.5  110.5  18  15  10  43
   88.3  92.8  96.8  13  11  8  10
   96.0  103.3  100.5  7  15  11  27
   94.3  94.5  99.0  10  12  11  19
   106.5  121.5  110.5  18  17  10  42
   106.5  115.5  107.0  8  13  14  47
   92.0  99.5  103.5  18  16  8  18
   102.0  99.8  103.3  13  12  14  28
   108.3  122.3  108.5  15  19  12  41
   106.8  119.0  106.8  14  20  12  37
   102.5  109.3  103.8  9  17  13  32
   92.5  102.5  99.3  13  15  6  23
   102.8  113.8  106.8  17  20  10  32
   83.3  87.3  96.3  1  5  9  15
   94.8  101.8  99.8  7  16  11  24
   103.5  112.0  110.8  18  13  12  37
   89.5  96.0  97.3  7  15  11  14
   84.3  89.8  94.3  8  8  8  9
   104.3  109.5  106.5  14  12  12  36
   106.0  118.5  105.0  12  16  11  39];  

   z = zx(:, 1:3);
   x = zx(:,4:7);
   [Bs, Cs, BL, CL]= SCA(x,z,3)