function printinfo(opt)  fprintf(opt.st,'Sammon mapping, error function with the parameter q=%d\n',opt.q);  switch (opt.optim)    case 'pn',        fprintf(opt.st,'Minimization by Pseudo-Newton algorithm\n');    case 'scg',        fprintf(opt.st,'Minimization by Scaled Conjugate Gradients algorithm\n');    otherwise       error(strcat('Possible initialization methods: pn (Pseudo-Newton), ',...                   'or scg (Scaled Conjugate Gradients).'));  endreturn% **********************************************************************************% [I,J,P] = MDS_REPPOINTS(D)%% Finds the indices of repeated/left points. J contains the indices of% repeated points in D. This means that for each j in J, there is a point% k ~= j such that D(J(j),k) = 0. I contains the indices of the remaining % points, and P those of the points in I that were identical to those in J.% Directly used in MDS routine.function [I,J,P] = mds_reppoints(D)  epsilon = 1e-20;      % Differences smaller than this are assumed to be zero.  [m,mm] = size(D);  I = 1:m; J = []; P = [];  if (m == mm) & (all(abs(D(1:m+1:end)) <= epsilon))    K = intersect (find (triu(ones(m),1)), find(D < epsilon));    if (~isempty(K))      P = mod(K,m);      J = fix(K./m) + 1;                 I(J) = [];                        end  else    [J,P] = find(D<=epsilon);     I(J)  = [];  endreturn;% **********************************************************************************% MDS_SETOPT  Sets parameters for the MDS mapping%%   OPT = MDS_SETOPT(OPT_GIVEN)% % INPUT%   OPT_GIVEN Parameters for the MDS mapping, described below (default: [])%% OUTPUT%   OPT       Structure of chosen options; if OPT_GIVEN is empty, then %               OPT.q       = 0%               OPT.optim   = 'pn'%               OPT.init    = 'cs'%               OPT.etol    = 1e-6 (the precise value depends on q)%               OPT.maxiter = inf %               OPT.itmap   = 'yes' %               OPT.isratio = 0 %               OPT.st      = 1 %               OPT.inspect = 2 %% DESCRIPTION  % Parameters for the MDS mapping can be set or changed. OPT_GIVEN consists% of the following fields: 'q', 'init', 'optim', 'etol', 'maxiter','itmap',% 'isratio', 'st' and 'inspect'. OPTIONS can include all the fields or some% of them only. The fields of OPT have some default values, which can be% changed by the OPT_GIVEN field values. If OPT_GIVEN is empty, then OPT% contains all default values. For a description of the fields, see MDS.%% Copyright: Elzbieta Pekalska, ela@ph.tn.tudelft.nl, 2000-2003% Faculty of Applied Sciences, Delft University of Technology%function opt = mds_setopt (opt_given)  opt.q       = 0;  opt.init    = 'cs';   opt.optim   = 'pn';  opt.st      = 1;   opt.itmap   = 'yes';   opt.maxiter = inf;   opt.inspect = 2;   opt.etol    = inf;   opt.isratio = 0;   opt.nanid   = [];    % Here are some extra values; set up in the MDS routine.  opt.plotlab = [];    % Not to be changed by the user.  if (~isempty(opt_given))    if (~isstruct(opt_given))      error('OPTIONS should be a structure with at least one of the following fields: q, init, etol, optim, maxiter, itmap, isratio, st or inspect.');    end    fn = fieldnames(opt_given);    if (~all(ismember(fn,fieldnames(opt))))      error('Wrong field names; valid field names are: q, init, optim, etol, maxiter, itmap, isratio, st or inspect.')    end     for i = 1:length(fn)      opt = setfield(opt,fn{i},getfield(opt_given,fn{i}));    end  end  if (isempty(intersect(opt.q,-2:1:2)))    error ('OPTIONS.q should be -2, -1, 0, 1 or 2.');  end  if (opt.maxiter < 2)    error ('OPTIONS.iter should be at least 1.');  end  if (isinf(opt.etol))    switch (opt.q)                 % Different defaults for different stresses.      case -2,         opt.etol = 1e-6;      case {-1,0}        opt.etol = 10*sqrt(eps);       case {1,2}        opt.etol = sqrt(eps);       end   elseif (opt.etol <= 0) | (opt.etol >= 0.1)    error ('OPTIONS.etol should be positive and smaller than 0.1.');  end  if (~ismember(opt.optim, {'pn','scg'}))    error('OPTIONS.optim should be pn or scg.');  end  if (~ismember(opt.itmap, {'yes','no'}))    error('OPTIONS.itmap should be yes or no.');  endreturn% **********************************************************************************% MDS_SAMMAP Sammon iterative nonlinear mapping for MDS%%   [YMAP,ERR] = MDS_SAMMAP(D,Y,YREP,OPTIONS)% % INPUT%   D       Square (M x M) dissimilarity matrix%   Y       M x N matrix containing starting configuration, or%   YREP    Configuration of the representation points%   OPTIONS Various parameters of the  minimization procedure put into %           a structure consisting of the following fields: 'q', 'optim',%           'init','etol','maxiter', 'isratio', 'itmap', 'st' and 'inspect'%           (default:%             OPTIONS.q       = 0          %             OPTIONS.optim   = 'pn'      %             OPTIONS.init    = 'cs'%             OPTIONS.etol    = 1e-6 (the precise value depends on q)%             OPTIONS.maxiter = inf %             OPTIONS.isratio = 0 %             OPTIONS.itmap   = 'yes' %             OPTIONS.st      = 1 %             OPTIONS.inspect = 2).% % OUTPUT%   YMAP     Mapped configuration%   ERR     Sammon stress %% DESCRIPTION  % Maps the objects given by a symmetric distance matrix D (with a zero% diagonal) onto, say, an N-dimensional configuration YMAP by an iterative% minimization of a variant of the Sammon stress. The minimization starts% from the initial configuration Y; see MDS_INIT.%% YREP is the Sammon configuration of the representation set. It is used% when new points have to be projected. In other words, if D is an M x M% symmetric distance matrix, then YREP is empty; if D is an M x N matrix,% then YMAP is sought such that D can approximate the distances between YMAP% and YREP.%% Missing values can be handled by marking them by NaN in D.%% SEE ALSO% MAPPINGS, MDS, MDS_CS, MDS_INIT, MDS_SETOPT%% Copyright: Elzbieta Pekalska, ela@ph.tn.tudelft.nl, 2000-2003% Faculty of Applied Sciences, Delft University of Technology%function [y,err] = mds_sammap(Ds,y,yrep,opt)  if (nargin < 4)    opt = [];                 % Will be filled by MDS_SETOPT, below.  end  if isempty(opt), opt = mds_setopt(opt); end  if (nargin < 3)    yrep = [];   end  % Extract labels and calculate distance matrix.  [m,n] = size(y);  if (~isempty(yrep))    replab = getlab(yrep);    yrep = +yrep;    D = sqrt(distm(y,yrep));   else    D = sqrt(distm(y));         yrep = [];    replab = [];  end  if (isempty(opt.plotlab))    opt.plotlab = ones(m,1);  end  it   = 0;           % Iteration number.  eold = inf;         % Previous stress (error).  % Calculate initial stress.  [e,a]= mds_samstress(opt.q,y,yrep,Ds,D,opt.nanid,opt.isratio); err = e;  fprintf(opt.st,'iteration: %4i   stress: %3.8d\n',it,e);    switch (opt.optim)      case 'pn',        % Pseudo-Newton minimization.      epr   = e;      % Previous values of E, EOLD for line search algorithm.      eepr  = e;      add   = 1e-4;   % Avoid G2 equal 0; see below.            BETA   = 1e-3;   % Parameter for line search algorithm.      lam1   = 0;           lam2   = 1;      % LAM (the step factor) lies in (LAM1, LAM2).      LMIN   = 1e-10;  % Minimum acceptable value for LAM in line search.      lambest = 1e-4; % LAM = LAMBEST, if LAM was not found.      % Loop until the error change falls below the tolerance or until      % the maximum number of iterations is reached.      while (abs(e-eold) >= opt.etol*(1 + abs(e))) & (it < opt.maxiter)        % Plot progress if requested.        if (opt.st > 0) & (opt.inspect > 0) & (mod(it,opt.inspect) == 0)%          if (it == 0), figure(1); clf; end          mds_plot(y,yrep,opt.plotlab,replab,e);         end        yold = y; eold = e;        % Calculate the stress.        [e,a]      = mds_samstress (opt.q,yold,yrep,Ds,[],opt.nanid,opt.isratio);        % Calculate gradients and pseudo-Newton update P.        [g1,g2,cc] = mds_gradients (opt.q,yold,yrep,a*Ds,[],opt.nanid);        p        = (g1/cc)./(abs(g2/cc) + add);        slope    = g1(:)' * p(:);        lam = 1;                         % Search for a suitable step LAM using a line search.        while (1)          % Take a step and calculate the delta error DE.          y      = yold + lam .* p;            [e,a]  = mds_samstress(opt.q,y,yrep,Ds,[],opt.nanid,opt.isratio);          de     = e - eold;          % Stop if the condition for a suitable lam is fulfilled.          if (de < BETA * lam * slope)            break;          end             % Try to find a suitable step LAM.          if (lam ~= 1)            r1 = de - lam*slope;             r2 = epr - eepr - lam2*slope;             aa = (2*de + r1/lam^2 - r2/lam2^2)/(lam2-lam1)^3;            bb = ((-lam2*r1)/lam^2 +(lam*r2)/lam2^2)/(lam2-lam1)^2;            if (abs(aa) <= eps)              lamtmp = -0.5 * slope/bb;            else              lamtmp = (-bb + sqrt(max(0,(bb^2 - 3*aa*slope))))/(3*aa);              end            % Prevent LAM from becoming too large.            if (lamtmp > 0.5 * lam), lamtmp = 0.5 * lam; end          else            lamtmp = -0.5 * slope/(e - eold - slope);          end          % Prevent LAM from becoming too small.          lam2 = lam;              lam  = max(lamtmp, 0.1*lam);            epr  = e;          eepr = eold;             if (lam < LMIN)            y = yold + lambest .* p;            [e,a] = mds_samstress(opt.q,y,yrep,Ds,[],opt.nanid,opt.isratio);            break;          end        end        it  = it + 1;         err = [err; e];        fprintf(opt.st,'iteration: %4i   stress: %3.8d\n',it,e);       end    case 'scg',              % Scaled Conjugate Gradient minimization.      sigma0  = 1e-8;                              lambda  = 1e-8;       % Regularization parameter.      lambda1 = 0;                                                                   % Calculate initial stress and direction of decrease P.      [e,a]   = mds_samstress(opt.q,y,yrep,Ds,D,opt.nanid,opt.isratio);         g1      = mds_gradients(opt.q,y,yrep,a*Ds,D,opt.nanid);   % Gradient      p       = -g1;                                            % Direction of decrease      gnorm2  = g1(:)' * g1(:);      success = 1;                                                    % Sanity check.      if (gnorm2 < 1e-15)        prwarning(2,['Gradient is nearly zero: ' gnorm2 ' (unlikely); initial configuration is the right one.']);          return      end      % Loop until the error change falls below the tolerance or until      % the maximum number of iterations is reached.      while (abs(eold - e) > opt.etol * (1 + e)) & (it < opt.maxiter)        g0     = g1;            % Previous gradient        pnorm2 = p(:)' * p(:);        eold   = e;