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