% Plot progress if requested.        if (opt.inspect > 0) & (mod(it,opt.inspect) == 0) %          if (it == 0), figure(1); clf; end          mds_plot(y,yrep,opt.plotlab,replab,e);         end        if (success)          sigma = sigma0/sqrt(pnorm2);      % SIGMA: a small step from y to yy          yy    = y + sigma .*p;          [e,a] = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio);           g2    = mds_gradients(opt.q,yy,yrep,a*Ds,[],opt.nanid); % Gradient for yy          s     = (g2-g1)/sigma;        % Approximation of Hessian*P directly, instead of computing          delta = p(:)' * s(:);         % the Hessian, since only DELTA = P'*Hessian*P is in fact needed        end        % Regularize the Hessian indirectly to make it positive definite.        % DELTA is now computed as P'* (Hessian + regularization * Identity) * P        delta = delta + (lambda1 - lambda) * pnorm2;          % Indicate if the Hessian is negative definite; the regularization above was too small.        if (delta < 0)          lambda1 = 2 * (lambda - delta/pnorm2);  % Now the Hessian will be positive definite          delta   = -delta + lambda * pnorm2;     % This is obtained after plugging lambda1           lambda  = lambda1;                      % into the formulation a few lines above.        end        mi = - p(:)' * g1(:);        yy = y + (mi/delta) .*p;                  % mi/delta is a step size         [ee,a] = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio);         % This minimization procedure is based on the second order approximation of the          % stress function by using the gradient and the Hessian approximation. The Hessian         % is regularized, but maybe not sufficiently. The ratio Dc (<=1) below indicates         % a proper approximation, if Dc is close to 1.        Dc = 2 * delta/mi^2 * (e - ee);         e  = ee;        % If Dc > 0, then the stress can be successfully decreased.        success = (Dc >= 0);        if (success)          y       = yy;          [ee,a]  = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio);           g1      = mds_gradients(opt.q,yy,yrep,a*Ds,[],opt.nanid);           gnorm2  = g1(:)' * g1(:);          lambda1 = 0;          beta = max(0,(g1(:)'*(g1(:)-g0(:)))/mi);            p   = -g1 + beta .* p;                 % P is a new conjugate direction            if (g1(:)'*p(:) >= 0 | mod(it-1,n*m) == 0),            p = -g1;                            % No much improvement, restart          end             if (Dc >= 0.75)            lambda = 0.25 * lambda;             % Make the regularization smaller          end          it  = it + 1;           err = [err; e];          fprintf (opt.st,'iteration: %4i   stress: %3.8d\n',it,e);         else        % Dc < 0            % Note that for Dc < 0, the iteration number IT is not increased,           % so the Hessian is further regularized until SUCCESS equals 1.          lambda1 = lambda;        end      % The approximation of the Hessian was poor or the stress was not       % decreased (Dc < 0), hence the regularization lambda is enlarged.      if (Dc < 0.25)        lambda = lambda + delta * (1 - Dc)/pnorm2;      end    end  endreturn;% **********************************************************************************%MDS_STRESS - Calculate Sammon stress during optimization%%   E = MDS_SAMSTRESS(Q,Y,YREP,Ds,D,NANINDEX)%% INPUT%   Q          Indicator of the Sammon stress; Q = -2,-1,0,1,2%   Y         Current lower-dimensional configuration%   YREP      Configuration of the representation objects; it should be%             empty when no representation set is considered%   Ds        Original distance matrix%   D         Approximate distance matrix (optional; otherwise computed from Y)%   NANINDEX  Index of the missing values; marked in Ds by NaN (optional; to%               be found in Ds)%% OUTPUT%   E         Sammon stress%% DESCRIPTION% Computes the Sammon stress between the original distance matrix Ds and the% approximated distance matrix D between the mapped configuration Y and the% configuration of the representation set YREP, expressed as follows:%% E = 1/(sum_{i<j} Ds_{ij}^(q+2)) sum_{i<j} (Ds_{ij} - D_{ij})^2 * Ds_{ij}^q%% It is directly used in the MDS_SAMMAP routine.%% Copyright: Elzbieta Pekalska, Robert P.W. Duin, ela@ph.tn.tudelft.nl, 2000-2003% Faculty of Applied Sciences, Delft University of Technology%function [e,alpha] = mds_samstress (q,y,yrep,Ds,D,nanindex,isratio)  % If D is not given, calculate it from Y, the current mapped points.  if (nargin < 5) | (isempty(D))    if (~isempty(yrep))      D = sqrt(distm(y,yrep));         else      D = sqrt(distm(y));         end  end  if (nargin < 6)    nanindex = [];    % Not given, so calculate below.  end  if (nargin < 7)    isratio = 0;      % Assume this is meant.  end  todefine = isempty(yrep);  [m,n]  = size(y); [mm,k] = size(Ds);  if (m ~= mm)    error('The sizes of Y and Ds do not match.');  end  if (any(size(D) ~= size(Ds)))    error ('The sizes of D and Ds do not match.');  end  m2 = m*k;  % Convert to double.  D  = +D; Ds = +Ds;  % I is the index of non-NaN, non-zero (> eps) values to be included   % for the computation of the stress.  I = 1:m2;   if (~isempty(nanindex))    I(nanindex) = [];  end  O = [];  if (todefine), O = 1:m+1:m2; end                                 I = setdiff(I,O);  % If OPTIONS.isratio is set, calculate optimal ALPHA to scale with.  if (isratio)    alpha = sum((Ds(I).^q).*D(I).^2)/sum((Ds(I).^(q+1)).*D(I));    Ds    = alpha*Ds;  else    alpha = 1;   end      % C is the normalization factor.  c = sum(Ds(I).^(q+2));   % If Q == 0, prevent unnecessary calculation (X^0 == 1).  if (q ~= 0)    e = sum(Ds(I).^q .*((Ds(I)-D(I)).^2))/c;  else    e = sum(((Ds(I)-D(I)).^2))/c;  endreturn; % **********************************************************************************%MDS_GRADIENTS - Gradients for variants of the Sammon stress%%   [G1,G2,CC] = MDS_GRADIENTS(Q,Y,YREP,Ds,D,NANINDEX)%% INPUT%   Q          Indicator of the Sammon stress; Q = -2,-1,0,1,2%   Y         Current lower-dimensional configuration%   YREP      Configuration of the representation objects; it should be%             empty when no representation set is considered%   Ds        Original distance matrix%   D         Approximate distance matrix (optional; otherwise computed from Y)%   nanindex  Index of missing values; marked in Ds by NaN (optional;%             otherwise found from Ds)%% OUTPUT%   G1        Gradient direction%   G2        Approximation of the Hessian by its diagonal  %% DESCRIPTION  % This is a routine used directly in the MDS_SAMMAP routine.%% SEE ALSO% MDS, MDS_INIT, MDS_SAMMAP, MDS_SAMSTRESS%% Copyright: Elzbieta Pekalska, Robert P.W. Duin, ela@ph.tn.tudelft.nl, 2000-2003% Faculty of Applied Sciences, Delft University of Technology%function [g1,g2,c] = mds_gradients(q,y,yrep,Ds,D,nanindex)  % If D is not given, calculate it from Y, the current mapped points.  if (nargin < 5) | (isempty(D))    if (~isempty(yrep))      D = sqrt(distm(y,yrep));         else      D = sqrt(distm(y));         end  end  % If NANINDEX is not given, find it from Ds.  if (nargin < 6)    nanindex = find(isnan(Ds(:))==1);    end  % YREP is empty if no representation set is defined yet.  % This happens when the mapping should be defined.  todefine = (isempty(yrep));  if (todefine)     yrep  = y;  end  [m,n]  = size(y); [mm,k] = size(Ds);  if (m ~= mm)    error('The sizes of Y and Ds do not match.');  end  if (any(size(D) ~= size(Ds)))    error ('The sizes of D and Ds do not match.');  end  m2 = m*k;  % Convert to doubles.  D  = +D; Ds = +Ds;  % I is the index of non-NaN, non-zero (> eps) values to be included   % for the computation of the gradient and the Hessians diagonal.  I = (1:m2)';  if (~isempty(nanindex))    I(nanindex) = [];  end  K    = find(Ds(I) <= eps | D(I) <= eps);  I(K) = [];  % C is a normalization factor.  c = -2/sum(Ds(I).^(q+2));   % Compute G1, the gradient.  h1 = zeros(m2,1);  if (q == 0)                    % Prevent unnecessary computation when Q = 0.    h1(I) = (Ds(I)-D(I)) ./ D(I);    else      h1(I) = (Ds(I)-D(I)) ./ (D(I).*(Ds(I).^(-q)));  end  h1 = reshape (h1',m,k);  g2 = h1 * ones(k,n);          % Here G2 is assigned only temporarily,  g1 = c * (g2.*y - h1*yrep);    % for the computation of G1.       % Compute G2, the diagonal of the Hessian, if requested.  if (nargout > 1)    h2 = zeros(m2,1);    switch (q)      case -2,         h2(I) = -1./(Ds(I).*D(I).^3);      case -1,         h2(I) = -1./(D(I).^3);      case 0,         h2(I) = - Ds(I)./(D(I).^3);      case 1,         h2(I) = - Ds(I).^2./(D(I).^3);      case 2,         h2(I) = -(Ds(I)./D(I)).^3;      end     h2 = reshape (h2',m,k);    g2 = c * (g2 + (h2*ones(k,n)).*y.^2 + h2*yrep.^2 - 2*(h2*yrep).*y);    endreturn;% **********************************************************************************%MDS_PLOT Plots the results of the MDS mapping in a 2D or 3D figure% %   MDS_PLOT(Y,YREP,LAB,REPLAB,E)% % INPUT%   Y       Configuration in 2D or 3D space%   YREP    Configuration of the representation points in 2D or 3D space%   LAB     Labels of Y%   REPLAB  Labels of YREP%   E       Stress value%% DESCRIPTION  % Used directly in MDS_SAMMAP routine.function mds_plot (y,yrep,lab,replab,e)  % This is done in a rather complicated way in order to speed up drawing.  K = min(size(y,2),3);  if (K > 1)    y   = +y;    col = 'brmk';    sym = ['+*xosdv^<>p']';    ii  = [1:44];    s   = [col(ii-floor((ii-1)/4)*4)' sym(ii-floor((ii-1)/11)*11)];    [nlab,lablist] = renumlab([replab;lab]); c = max(nlab);    [ms,ns] = size(s); if (ms == 1), s = setstr(ones(m,1)*s); end    yy = [yrep; y];      cla;     if (K == 2)      hold on;      for i = 1:c        J = find(nlab==i); plot(yy(J,1),yy(J,2),s(i,:));      end    else       for i = 1:c        J = find(nlab==i); plot3(yy(J,1),yy(J,2),yy(J,3),s(i,:));        if (i == 1), hold on; grid on; end      end    end    title(sprintf('STRESS: %f', e));    axis equal;      drawnow;  endreturn;