function [Y, R, E] = Isomap(D, n_fcn, n_size, options); % ISOMAP   Computes Isomap embedding using the algorithm of %             Tenenbaum, de Silva, and Langford (2000). %% [Y, R, E] = isomap(D, n_fcn, n_size, options); %% Input:%    D = N x N matrix of distances (where N is the number of data points)%    n_fcn = neighborhood function ('epsilon' or 'k') %    n_size = neighborhood size (value for epsilon or k) %%    options.dims = (row) vector of embedding dimensionalities to use%                        (1:10 = default)%    options.comp = which connected component to embed, if more than one. %                        (1 = largest (default), 2 = second largest, ...)%    options.display = plot residual variance and 2-D embedding?%                        (1 = yes (default), 0 = no)%    options.overlay = overlay graph on 2-D embedding?  %                        (1 = yes (default), 0 = no)%    options.verbose = display progress reports? %                        (1 = yes (default), 0 = no)%% Output: %    Y = Y.coords is a cell array, with coordinates for d-dimensional embeddings%         in Y.coordsdx33zjb.  Y.index contains the indices of the points embedded.%    R = residual variances for embeddings in Y%    E = edge matrix for neighborhood graph%%    BEGIN COPYRIGHT NOTICE%%    Isomap code -- (c) 1998-2000 Josh Tenenbaum%%    This code is provided as is, with no guarantees except that %    bugs are almost surely present.  Published reports of research %    using this code (or a modified version) should cite the %    article that describes the algorithm: %%      J. B. Tenenbaum, V. de Silva, J. C. Langford (2000).  A global%      geometric framework for nonlinear dimensionality reduction.  %      Science 290 (5500): 2319-2323, 22 December 2000.  %%    Comments and bug reports are welcome.  Email to jbt@psych.stanford.edu. %    I would also appreciate hearing about how you used this code, %    improvements that you have made to it, or translations into other%    languages.    %%    You are free to modify, extend or distribute this code, as long %    as this copyright notice is included whole and unchanged.  %%    END COPYRIGHT NOTICE%%%%% Step 0: Initialization and Parameters %%%%%N = size(D,1); if ~(N==size(D,2))     error('D must be a square matrix'); end; if n_fcn=='k'     K = n_size;      if ~(K==round(K))         error('Number of neighbors for k method must be an integer');     endelseif n_fcn=='epsilon'     epsilon = n_size; else      error('Neighborhood function must be either epsilon or k'); endif nargin < 3     error('Too few input arguments'); elseif nargin < 4     options = struct('dims',1:10,'overlay',1,'comp',1,'display',1,'verbose',1); endinf =  1000*max(max(D))*N;  %% effectively infinite distanceif ~isfield(options,'dims')     options.dims = 1:10; endif ~isfield(options,'overlay')     options.overlay = 1; endif ~isfield(options,'comp')     options.comp = 1; endif ~isfield(options,'display')     options.display = 1; endif ~isfield(options,'verbose')     options.verbose = 1; enddims = options.dims; comp = options.comp; overlay = options.overlay; displ = options.display; verbose = options.verbose; Y.coords = cell(length(dims),1); R = zeros(1,length(dims)); %%%%% Step 1: Construct neighborhood graph %%%%%disp('Constructing neighborhood graph...'); if n_fcn == 'k'     [tmp, ind] = sort(D);      for i=1:N          D(i,ind((2+K):end,i)) = inf;      endelseif n_fcn == 'epsilon'     warning off    %% Next line causes an unnecessary warning, so turn it off     D =  D./(D<=epsilon);      D = min(D,inf);      warning onendD = min(D,D');    %% Make sure distance matrix is symmetricif (overlay == 1)     E = int8(1-(D==inf));  %%  Edge information for subsequent graph overlayend% Finite entries in D now correspond to distances between neighboring points. % Infinite entries (really, equal to INF) in D now correspond to %   non-neighoring points. %%%%% Step 2: Compute shortest paths %%%%%disp('Computing shortest paths...'); % We use Floyd's algorithm, which produces the best performance in Matlab. % Dijkstra's algorithm is significantly more efficient for sparse graphs, % but requires for-loops that are very slow to run in Matlab.  A significantly % faster implementation of Isomap that calls a MEX file for Dijkstra's % algorithm can be found in isomap2.m (and the accompanying files% dijkstra.c and dijkstra.dll). tic; for k=1:N     D = min(D,repmat(D(:,k),[1 N])+repmat(D(k,:),[N 1]));      if ((verbose == 1) & (rem(k,20) == 0))           disp([' Iteration: ' num2str(k) 'Estimated time to completion: 'num2str((N-k)*toc/k/60) ' minutes']);      endend%%%%% Remove outliers from graph %%%%%disp('Checking for outliers...'); n_connect = sum(~(D==inf));        %% number of points each point connects to[tmp, firsts] = min(D==inf);       %% first point each point connects to[comps, I, J] = unique(firsts);    %% represent each connected component oncesize_comps = n_connect(comps);     %% size of each connected component[tmp, comp_order] = sort(size_comps);  %% sort connected components by sizecomps = comps(comp_order(end:-1:1));    size_comps = size_comps(comp_order(end:-1:1)); n_comps = length(comps);               %% number of connected componentsif (comp>n_comps)                     comp=1;                              %% default: use largest componentenddisp(['Number of connected components in graph: ' num2str(n_comps)]); disp(['Embedding component ' num2str(comp) ' with ' num2str(size_comps(comp)) ' points.']); Y.index = find(firsts==comps(comp)); D = D(Y.index, Y.index); N = length(Y.index); %%%%% Step 3: Construct low-dimensional embeddings (Classical MDS) %%%%%disp('Constructing low-dimensional embeddings (Classical MDS)...'); opt.disp = 0; [vec, val] = eigs(-.5*(D.^2 - sum(D.^2)'*ones(1,N)/N - ones(N,1)*sum(D.^2)/N + sum(sum(D.^2))/(N^2)), max(dims), 'LR', opt); h = real(diag(val)); [foo,sorth] = sort(h);  sorth = sorth(end:-1:1); val = real(diag(val(sorth,sorth))); vec = vec(:,sorth); D = reshape(D,N^2,1); for di = 1:length(dims)     if (dims(di)<=N)         Y.coords{di} = real(vec(:,1:dims(di)).*(ones(N,1)*sqrt(val(1:dims(di)))'))';          r2 = 1-corrcoef(reshape(real(L2_distance(Y.coords{di}, Y.coords{di})),N^2,1),D).^2;          R(di) = r2(2,1);          if (verbose == 1)             disp(['Isomap on 'num2str(N) 'points with dimensionality ' num2str(dims(di)) '  --> residual variance = ' num2str(R(di))]);          end     endendclear D; %%%%%%%%%%%%%%%%%% Graphics %%%%%%%%%%%%%%%%%%if (displ==1)     %%%%% Plot fall-off of residual variance with dimensionality %%%%%     figure;     hold on     plot(dims, R, 'bo');      plot(dims, R, 'b-');      hold off     ylabel('Residual variance');      xlabel('Isomap dimensionality');      %%%%% Plot two-dimensional configuration %%%%%     twod = find(dims==2);      if ~isempty(twod)         figure;         hold on;         plot(Y.coords{twod}(1,:), Y.coords{twod}(2,:), 'ro');          if (overlay == 1)             gplot(E(Y.index, Y.index), [Y.coords{twod}(1,:); Y.coords{twod}(2,:)]');              title('Two-dimensional Isomap embedding (with neighborhood graph).');          else             title('Two-dimensional Isomap.');          end         hold off;     endendreturn;