wolffd@0: function [P] = cca(D, P, epochs, Mdist, alpha0, lambda0)
wolffd@0: 
wolffd@0: %CCA Projects data vectors using Curvilinear Component Analysis.
wolffd@0: %
wolffd@0: % P = cca(D, P, epochs, [Dist], [alpha0], [lambda0])
wolffd@0: %
wolffd@0: %  P = cca(D,2,10);           % projects the given data to a plane
wolffd@0: %  P = cca(D,pcaproj(D,2),5); % same, but with PCA initialization
wolffd@0: %  P = cca(D, 2, 10, Dist);   % same, but the given distance matrix is used
wolffd@0: %  
wolffd@0: %  Input and output arguments ([]'s are optional):
wolffd@0: %   D          (matrix) the data matrix, size dlen x dim
wolffd@0: %              (struct) data or map struct            
wolffd@0: %   P          (scalar) output dimension
wolffd@0: %              (matrix) size dlen x odim, the initial projection
wolffd@0: %   epochs     (scalar) training length
wolffd@0: %   [Dist]     (matrix) pairwise distance matrix, size dlen x dlen.
wolffd@0: %                       If the distances in the input space should
wolffd@0: %                       be calculated otherwise than as euclidian
wolffd@0: %                       distances, the distance from each vector
wolffd@0: %                       to each other vector can be given here,
wolffd@0: %                       size dlen x dlen. For example PDIST
wolffd@0: %                       function can be used to calculate the
wolffd@0: %                       distances: Dist = squareform(pdist(D,'mahal'));
wolffd@0: %   [alpha0]   (scalar) initial step size, 0.5 by default
wolffd@0: %   [lambda0]  (scalar) initial radius of influence, 3*max(std(D)) by default
wolffd@0: %  
wolffd@0: %   P          (matrix) size dlen x odim, the projections
wolffd@0: %
wolffd@0: % Unknown values (NaN's) in the data: projections of vectors with
wolffd@0: % unknown components tend to drift towards the center of the
wolffd@0: % projection distribution. Projections of totally unknown vectors are
wolffd@0: % set to unknown (NaN).
wolffd@0: %
wolffd@0: % See also SAMMON, PCAPROJ. 
wolffd@0: 
wolffd@0: % Reference: Demartines, P., Herault, J., "Curvilinear Component
wolffd@0: %   Analysis: a Self-Organizing Neural Network for Nonlinear
wolffd@0: %   Mapping of Data Sets", IEEE Transactions on Neural Networks, 
wolffd@0: %   vol 8, no 1, 1997, pp. 148-154.
wolffd@0: 
wolffd@0: % Contributed to SOM Toolbox 2.0, February 2nd, 2000 by Juha Vesanto
wolffd@0: % Copyright (c) by Juha Vesanto
wolffd@0: % http://www.cis.hut.fi/projects/somtoolbox/
wolffd@0: 
wolffd@0: % juuso 171297 040100
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% Check arguments 
wolffd@0: 
wolffd@0: error(nargchk(3, 6, nargin)); % check the number of input arguments
wolffd@0: 
wolffd@0: % input data
wolffd@0: if isstruct(D), 
wolffd@0:   if strcmp(D.type,'som_map'), D = D.codebook; else D = D.data; end
wolffd@0: end
wolffd@0: [noc dim] = size(D);
wolffd@0: noc_x_1  = ones(noc, 1); % used frequently
wolffd@0: me = zeros(1,dim); st = zeros(1,dim);
wolffd@0: for i=1:dim,
wolffd@0:   me(i) = mean(D(find(isfinite(D(:,i))),i));
wolffd@0:   st(i) = std(D(find(isfinite(D(:,i))),i));
wolffd@0: end
wolffd@0: 
wolffd@0: % initial projection
wolffd@0: if prod(size(P))==1, 
wolffd@0:   P = (2*rand(noc,P)-1).*st(noc_x_1,1:P) + me(noc_x_1,1:P); 
wolffd@0: else
wolffd@0:   % replace unknown projections with known values
wolffd@0:   inds = find(isnan(P)); P(inds) = rand(size(inds));
wolffd@0: end
wolffd@0: [dummy odim] = size(P);
wolffd@0: odim_x_1  = ones(odim, 1); % this is used frequently
wolffd@0: 
wolffd@0: % training length
wolffd@0: train_len = epochs*noc;
wolffd@0: 
wolffd@0: % random sample order
wolffd@0: rand('state',sum(100*clock));
wolffd@0: sample_inds = ceil(noc*rand(train_len,1));
wolffd@0: 
wolffd@0: % mutual distances
wolffd@0: if nargin<4 | isempty(Mdist) | all(isnan(Mdist(:))),
wolffd@0:   fprintf(2, 'computing mutual distances\r');
wolffd@0:   dim_x_1 = ones(dim,1);
wolffd@0:   for i = 1:noc,
wolffd@0:     x = D(i,:); 
wolffd@0:     Diff = D - x(noc_x_1,:);
wolffd@0:     N = isnan(Diff);
wolffd@0:     Diff(find(N)) = 0; 
wolffd@0:     Mdist(:,i) = sqrt((Diff.^2)*dim_x_1);
wolffd@0:     N = find(sum(N')==dim); %mutual distance unknown
wolffd@0:     if ~isempty(N), Mdist(N,i) = NaN; end
wolffd@0:   end
wolffd@0: else
wolffd@0:   % if the distance matrix is output from PDIST function
wolffd@0:   if size(Mdist,1)==1, Mdist = squareform(Mdist); end
wolffd@0:   if size(Mdist,1)~=noc, 
wolffd@0:     error('Mutual distance matrix size and data set size do not match'); 
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: % alpha and lambda
wolffd@0: if nargin<5 | isempty(alpha0) | isnan(alpha0), alpha0 = 0.5; end
wolffd@0: alpha = potency_curve(alpha0,alpha0/100,train_len);
wolffd@0: 
wolffd@0: if nargin<6 | isempty(lambda0) | isnan(lambda0), lambda0 = max(st)*3; end
wolffd@0: lambda = potency_curve(lambda0,0.01,train_len);
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% Action
wolffd@0: 
wolffd@0: k=0; fprintf(2, 'iterating: %d / %d epochs\r',k,epochs);
wolffd@0: 
wolffd@0: for i=1:train_len, 
wolffd@0:   
wolffd@0:   ind = sample_inds(i);     % sample index
wolffd@0:   dx = Mdist(:,ind);        % mutual distances in input space
wolffd@0:   known = find(~isnan(dx)); % known distances
wolffd@0: 
wolffd@0:   if ~isempty(known),
wolffd@0:     % sample vector's projection
wolffd@0:     y = P(ind,:);                 
wolffd@0: 
wolffd@0:     % distances in output space
wolffd@0:     Dy = P(known,:) - y(noc_x_1(known),:); 
wolffd@0:     dy = sqrt((Dy.^2)*odim_x_1);           
wolffd@0:   
wolffd@0:     % relative effect
wolffd@0:     dy(find(dy==0)) = 1;        % to get rid of div-by-zero's
wolffd@0:     fy = exp(-dy/lambda(i)) .* (dx(known) ./ dy - 1);
wolffd@0: 
wolffd@0:     % Note that the function F here is e^(-dy/lambda)) 
wolffd@0:     % instead of the bubble function 1(lambda-dy) used in the 
wolffd@0:     % paper.
wolffd@0:     
wolffd@0:     % Note that here a simplification has been made: the derivatives of the
wolffd@0:     % F function have been ignored in calculating the gradient of error
wolffd@0:     % function w.r.t. to changes in dy.
wolffd@0:     
wolffd@0:     % update
wolffd@0:     P(known,:) = P(known,:) + alpha(i)*fy(:,odim_x_1).*Dy;
wolffd@0:   end
wolffd@0: 
wolffd@0:   % track
wolffd@0:   if rem(i,noc)==0, 
wolffd@0:     k=k+1; fprintf(2, 'iterating: %d / %d epochs\r',k,epochs);
wolffd@0:   end
wolffd@0: 
wolffd@0: end
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% clear up
wolffd@0: 
wolffd@0: % calculate error
wolffd@0: error = cca_error(P,Mdist,lambda(train_len));
wolffd@0: fprintf(2,'%d iterations, error %f          \n', epochs, error);
wolffd@0: 
wolffd@0: % set projections of totally unknown vectors as unknown
wolffd@0: unknown = find(sum(isnan(D)')==dim);
wolffd@0: P(unknown,:) = NaN;
wolffd@0: 
wolffd@0: return;
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% tips
wolffd@0: 
wolffd@0: % to plot the results, use the code below
wolffd@0: 
wolffd@0: %subplot(2,1,1), 
wolffd@0: %switch(odim), 
wolffd@0: %  case 1, plot(P(:,1),ones(dlen,1),'x')
wolffd@0: %  case 2, plot(P(:,1),P(:,2),'x'); 
wolffd@0: %  otherwise, plot3(P(:,1),P(:,2),P(:,3),'x'); rotate3d on
wolffd@0: %end
wolffd@0: %subplot(2,1,2), dydxplot(P,Mdist);
wolffd@0: 
wolffd@0: % to a project a new point x in the input space to the output space
wolffd@0: % do the following:
wolffd@0: 
wolffd@0: % Diff = D - x(noc_x_1,:); Diff(find(isnan(Diff))) = 0; 
wolffd@0: % dx = sqrt((Diff.^2)*dim_x_1);
wolffd@0: % p = project_point(P,x,dx); % this function can be found from below
wolffd@0: % tlen = size(p,1);
wolffd@0: % plot(P(:,1),P(:,2),'bx',p(tlen,1),p(tlen,2),'ro',p(:,1),p(:,2),'r-')
wolffd@0: 
wolffd@0: % similar trick can be made to the other direction
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% subfunctions
wolffd@0: 
wolffd@0: function vals = potency_curve(v0,vn,l)
wolffd@0: 
wolffd@0:   % curve that decreases from v0 to vn with a rate that is 
wolffd@0:   % somewhere between linear and 1/t
wolffd@0:   vals = v0 * (vn/v0).^([0:(l-1)]/(l-1));
wolffd@0: 
wolffd@0: 
wolffd@0: function error = cca_error(P,Mdist,lambda)
wolffd@0: 
wolffd@0:   [noc odim] = size(P);
wolffd@0:   noc_x_1 = ones(noc,1);
wolffd@0:   odim_x_1 = ones(odim,1);
wolffd@0: 
wolffd@0:   error = 0;
wolffd@0:   for i=1:noc,
wolffd@0:     known = find(~isnan(Mdist(:,i)));
wolffd@0:     if ~isempty(known),   
wolffd@0:       y = P(i,:);                 
wolffd@0:       Dy = P(known,:) - y(noc_x_1(known),:);
wolffd@0:       dy = sqrt((Dy.^2)*odim_x_1);
wolffd@0:       fy = exp(-dy/lambda);
wolffd@0:       error = error + sum(((Mdist(known,i) - dy).^2).*fy);
wolffd@0:     end
wolffd@0:   end
wolffd@0:   error = error/2;
wolffd@0: 
wolffd@0: 
wolffd@0: function [] = dydxplot(P,Mdist)
wolffd@0: 
wolffd@0:   [noc odim] = size(P);
wolffd@0:   noc_x_1 = ones(noc,1);
wolffd@0:   odim_x_1 = ones(odim,1);
wolffd@0:   Pdist = zeros(noc,noc);
wolffd@0:     
wolffd@0:   for i=1:noc,
wolffd@0:     y = P(i,:);                 
wolffd@0:     Dy = P - y(noc_x_1,:);
wolffd@0:     Pdist(:,i) = sqrt((Dy.^2)*odim_x_1);
wolffd@0:   end
wolffd@0: 
wolffd@0:   Pdist = tril(Pdist,-1); 
wolffd@0:   inds = find(Pdist > 0); 
wolffd@0:   n = length(inds);
wolffd@0:   plot(Pdist(inds),Mdist(inds),'.');
wolffd@0:   xlabel('dy'), ylabel('dx')
wolffd@0: 
wolffd@0: 
wolffd@0: function p = project_point(P,x,dx)
wolffd@0: 
wolffd@0:   [noc odim] = size(P);
wolffd@0:   noc_x_1 = ones(noc,1);
wolffd@0:   odim_x_1 = ones(odim,1);
wolffd@0: 
wolffd@0:   % initial projection
wolffd@0:   [dummy,i] = min(dx);
wolffd@0:   y = P(i,:)+rand(1,odim)*norm(P(i,:))/20;
wolffd@0:  
wolffd@0:   % lambda 
wolffd@0:   lambda = norm(std(P));
wolffd@0: 
wolffd@0:   % termination
wolffd@0:   eps = 1e-3; i_max = noc*10;
wolffd@0:   
wolffd@0:   i=1; p(i,:) = y; 
wolffd@0:   ready = 0;
wolffd@0:   while ~ready,
wolffd@0: 
wolffd@0:     % mutual distances
wolffd@0:     Dy = P - y(noc_x_1,:);        % differences in output space
wolffd@0:     dy = sqrt((Dy.^2)*odim_x_1);  % distances in output space
wolffd@0:     f = exp(-dy/lambda);
wolffd@0:   
wolffd@0:     fprintf(2,'iteration %d, error %g \r',i,sum(((dx - dy).^2).*f));
wolffd@0: 
wolffd@0:     % all the other vectors push the projected one
wolffd@0:     fy = f .* (dx ./ dy - 1) / sum(f);
wolffd@0:   
wolffd@0:     % update    
wolffd@0:     step = - sum(fy(:,odim_x_1).*Dy);
wolffd@0:     y = y + step;
wolffd@0:   
wolffd@0:     i=i+1;
wolffd@0:     p(i,:) = y;   
wolffd@0:     ready = (norm(step)/norm(y) < eps | i > i_max);
wolffd@0: 
wolffd@0:   end
wolffd@0:   fprintf(2,'\n');
wolffd@0:      
wolffd@0: