D = D(Y.index, Y.index); 
N = length(Y.index); 
%%%%% Step 3: Construct 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},0)),N^2,1),D).^2; 
         R(di) = r2(2,1); 
     end
end
clear D; 


% --- L2_distance function
% Written by Roland Bunschoten, University of Amsterdam, 1999
function d = L2_distance(a,b,df)
if (size(a,1) == 1)
  a = [a; zeros(1,size(a,2))]; 
  b = [b; zeros(1,size(b,2))]; 
end
aa=sum(a.*a); bb=sum(b.*b); ab=a'*b; 
d = sqrt(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab);
d = real(d); 
if (df==1)
  d = d.*(1-eye(size(d)));
end


% --- HLLE function
% Written by David Donoho & Carrie Grimes, 2003.
function [Y, mse] = HLLE(X,k,d)
N = size(X,2);
if max(size(k)) ==1
    kvec = repmat(k,N,1);
elseif max(size(k)) == N
    kvec=k;
end;
%Compute Nearest neighbors
D1 = L2_distance(X,X,1);
dim = size(X,1);
nind = repmat(0, size(D1,1), size(D1,2));
dp = d*(d+1)/2;
W = repmat(0,dp*N,N);
if(mean(k)>d) 
  tol=1e-3; % regularlizer in case constrained fits are ill conditioned
else
  tol=0;
end;
for i=1:N
    tmp = D1(:,i);
    [ts, or] = sort(tmp);
    %take k nearest neighbors
    nind(or(2:kvec(i)+1),i) = 1;
    thisx = X(:,or(2:kvec(i)+1));
    %center using the mean 
    thisx = thisx - repmat(mean(thisx')',1,kvec(i));
    %compute local coordinates
    [U,D,Vpr] = svd(thisx);
    V = Vpr(:,1:d);
    %Neighborhood diagnostics
    vals = diag(D);
    mse(i) = sum(vals(d+1:end));
    %build Hessian estimator
    clear Yi; clear Pii;
    ct = 0;
    for mm=1:d
        startp = V(:,mm);
        for nn=1:length(mm:d)
            indles = mm:d;
            Yi(:,ct+nn) = startp.*(V(:,indles(nn)));
        end;
        ct = ct+length(mm:d);
    end;
    Yi = [repmat(1,kvec(i),1), V, Yi];
    %orthogonalize linear and quadratic forms
    [Yt, Orig] = mgs(Yi);
    Pii = Yt(:,d+2:end)';
    %double check weights sum to 1
    for j=1:dp
        if sum(Pii(j,:)) >0.0001
            tpp = Pii(j,:)./sum(Pii(j,:)); 
        else
            tpp = Pii(j,:);
        end;
        %fill weight matrix
       W((i-1)*dp+j, or(2:kvec(i)+1)) = tpp;
    end;
end;
%%%%%%%%%%%%%%%%%%%%Compute eigenanalysis of W
G=W'*W;
G = sparse(G);
options.disp = 0; 
options.isreal = 1; 
options.issym = 1;
tol=0;
[Yo,eigenvals] = eigs(G,d+1,tol,options);
Y = Yo(:,1:d)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0
%compute final coordinate alignment
R = Y'*Y;
R2 = R^(-1/2);
Y = Y*R2;


% --- leigs function for Laplacian eigenmap.
% Written by Belkin & Niyogi, 2002.
function [E,V] = leigs(DATA, TYPE, PARAM, NE) 
n = size(DATA,1);
A = sparse(n,n);
step = 100;  
for i1=1:step:n    
    i2 = i1+step-1;
    if (i2> n) 
      i2=n;
    end;
    XX= DATA(i1:i2,:);  
    dt = L2_distance(XX',DATA',0);
    [Z,I] = sort ( dt,2);
    for i=i1:i2
      for j=2:PARAM+1
	        A(i,I(i-i1+1,j))= Z(i-i1+1,j); 
	        A(I(i-i1+1,j),i)= Z(i-i1+1,j); 
      end;    
    end;
end;
W = A;
[A_i, A_j, A_v] = find(A);  % disassemble the sparse matrix
for i = 1: size(A_i)  
    W(A_i(i), A_j(i)) = 1;
end;
D = sum(W(:,:),2);   
L = spdiags(D,0,speye(size(W,1)))-W;
opts.tol = 1e-9;
opts.issym=1; 
opts.disp = 0; 
[E,V] = eigs(L,NE,'sm',opts);


% --- diffusionKernel function
% Written by R. Coifman & S. Lafon.
function [Y] = diffusionKernel (X,sigmaK,alpha,d)
D = L2_distance(X',X',1);
K = exp(-(D/sigmaK).^2);
p = sum(K);
p = p(:);
K1 = K./((p*p').^alpha);
v = sqrt(sum(K1));
v = v(:);
A = K1./(v*v');
if sigmaK >= 0.5
    thre = 1e-7;  
    M = max(max(A));
    A = sparse(A.*double(A>thre*M));
    [U,S,V] = svds(A,d+1);   %Sparse version.
    U = U./(U(:,1)*ones(1,d+1));
else
    [U,S,V] = svd(A,0);   %Full version.
    U = U./(U(:,1)*ones(1,size(U,1)));
end;
Y = U(:,2:d+1);


% --- mgs function: Modified Gram-Schmidt
% Used by HLLE function.
function [Q, R] = mgs(A);
[m, n] = size(A);   % Assume m>=n.
V = A;
R = zeros(n,n);
for i=1:n
    R(i,i) = norm(V(:,i));
    V(:,i) = V(:,i)/R(i,i);
    if (i < n)
        for j = i+1:n
            R(i,j) = V(:,i)' * V(:,j);
            V(:,j) = V(:,j) - R(i,j) * V(:,i);
        end;
     end;
 end;
 Q = V;


% --- function mdsFast for Multi-Dimensional Scaling
% Written by Michael D. Lee.
% Lee recommends metric=2, iterations=50, learnrate=0.05.
function [points]=mdsFast(d,dim)
[n, check] = size(d);
iterations = 30;
lr=0.05;   % learnrate
r=2;   % metric
% normalise distances to lie between 0 and 1
reshift=min(min(d));
d=d-reshift;
rescale=max(max(d));
d=d/rescale;
% calculate the variance of  the distance matrix
dbar=(sum(sum(d))-trace(d))/n/(n-1);
temp=(d-dbar*ones(n)).^2;
vard=.5*(sum(sum(temp))-trace(temp));
% initialize variables
its=0;
p=rand(n,dim)*.01-.005;
dh=zeros(n);
rinv=1/r;  %PT
kk=1:dim;  %PT 
% main loop for the given number of iterations
while (its<iterations)
   its=its+1;
   % randomly permute the objects to determine the order
   % in which they are pinned for this iteration
   pinning_order=randperm(n);
   for i=1:n
      m=pinning_order(i);
      % having pinned an object, move all of the other on each dimension
      % according to the learning rule   
      
      %PT: Vectorized the procedure, gives factor 6 speed up for n=100 dim=2
      indx=[1:m-1 m+1:n];                                                       
      pmat=repmat(p(m,:),[n 1])-p;                                              
      dhdum=sum(abs(pmat).^r,2).^rinv;
      dh(m,indx)=dhdum(indx)';
      dh(indx,m)=dhdum(indx);
      dhmat=lr*repmat((dhdum(indx)-d(m,indx)').*(dhdum(indx).^(1-r)),[1 dim]);
      p(indx,kk)=p(indx,kk)+dhmat.*abs(pmat(indx,kk)).^(r-1).*sign(pmat(indx,kk));
                    %plus sign in learning rule is due the sign of pmat
   end;
end;
points = p*rescale+reshift;


% --- Executes during object creation, after setting all properties.
function ExampleParamEdit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ExampleParamEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function ExampleParamEdit_Callback(hObject, eventdata, handles)
% hObject    handle to ExampleParamEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ExampleParamEdit as text
%        str2double(get(hObject,'String')) returns contents of ExampleParamEdit as a double


% --- Executes during object creation, after setting all properties.
function ParamEdit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ParamEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


% --- function pca
% PCA analysis of data set.
function [Y] = pca (X,d)
opts.disp = 0;
[Y,D] = eigs(X*X',d,'lm',opts);


function ParamEdit_Callback(hObject, eventdata, handles)


% --- Executes on button press in DiffKernelButton.
function DiffKernelButton_Callback(hObject, eventdata, handles)
handles.sigma = str2double(get(handles.SigmaEdit,'String'));
handles.sigma = max(0,abs(handles.sigma));
handles.alpha = str2double(get(handles.AlphaEdit,'String'));
handles.alpha = abs(handles.alpha);
handles.d = round(str2double(get(handles.DimEdit,'String')));
handles.d = max(1,handles.d);
updateString{1} = 'Running Diffusion Map.';
set(handles.UpdatesText,'String',updateString);
tic;
handles.Y = diffusionKernel(handles.X,handles.sigma,handles.alpha,handles.d);
runTime = toc;
guidata(hObject, handles);
PlotEmbedding(hObject,eventdata,handles,0);
assignin ('base','maniY',handles.Y);
updateString{2} = ['Diffusion Map complete: ',num2str(runTime),'s'];
updateString{3} = 'Embedding data written to matrix "maniY"';
set(handles.UpdatesText,'String',updateString);


% --- Executes during object creation, after setting all properties.
function AlphaEdit_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function AlphaEdit_Callback(hObject, eventdata, handles)



% --- Executes on button press in LTSAbutton.
function LTSAbutton_Callback(hObject, eventdata, handles)
handles.d = round(str2double(get(handles.DimEdit,'String')));
handles.d = max(1,handles.d);
handles.K = round(str2double(get(handles.KEdit,'String')));
handles.K = max(1,handles.K);
updateString{1} = 'Running LTSA.';
set(handles.UpdatesText,'String',updateString);
tic;
[T,NI] = LTSA(handles.X',handles.d,handles.K);
handles.Y = T';
runTime = toc;
guidata(hObject, handles);
PlotEmbedding(hObject,eventdata,handles,0);
assignin ('base','maniY',handles.Y);
updateString{2} = ['LTSA complete: ',num2str(runTime),'s'];
updateString{3} = 'Embedding data written to matrix "maniY"';
set(handles.UpdatesText,'String',updateString);
guidata(hObject, handles);


% --- LTSA function
% Written by Zhenyue Zhang & Hongyuan Zha, 2004.
% Reference: http://epubs.siam.org/sam-bin/dbq/article/41915
function [T,NI] = LTSA(data,d,K,NI)
[m,N] = size(data);  % m is the dimensionality of the input sample points. 
% Step 0:  Neighborhood Index
if nargin<4
    if length(K)==1
        K = repmat(K,[1,N]);
    end;
    NI = cell(1,N); 
    if m>N
        a = sum(data.*data); 
        dist2 = sqrt(repmat(a',[1 N]) + repmat(a,[N 1]) - 2*(data'*data));
        for i=1:N
            % Determine ki nearest neighbors of x_j
            [dist_sort,J] = sort(dist2(:,i));  
            Ii = J(1:K(i)); 
            NI{i} = Ii;
        end;
    else
        for i=1:N
            % Determine ki nearest neighbors of x_j
            x = data(:,i); ki = K(i);
            dist2 = sum((data-repmat(x,[1 N])).^2,1);    
            [dist_sort,J] = sort(dist2);  
            Ii = J(1:ki);  
            NI{i} = Ii;
        end;
    end;
else
    K = zeros(1,N);
    for i=1:N
        K(i) = length(NI{i});
    end;
end;
% Step 1:  local information
BI = {}; 
Thera = {}; 
for i=1:N
    % Compute the d largest right singular eigenvectors of the centered matrix
    Ii = NI{i}; ki = K(i);
    Xi = data(:,Ii)-repmat(mean(data(:,Ii),2),[1,ki]);
    W = Xi'*Xi; W = (W+W')/2;
    [Vi,Si] = schur(W);
    [s,Ji] = sort(-diag(Si)); 
    Vi = Vi(:,Ji(1:d));  
    % construct Gi
    Gi = [repmat(1/sqrt(ki),[ki,1]) Vi];  
    % compute the local orthogonal projection Bi = I-Gi*Gi' 
    % that has the null space span([e,Theta_i^T]). 
    BI{i} = eye(ki)-Gi*Gi';    
end;
B = speye(N);
for i=1:N
    Ii = NI{i};
    B(Ii,Ii) = B(Ii,Ii)+BI{i};
    B(i,i) = B(i,i)-1;
end;
B = (B+B')/2;
options.disp = 0; 
options.isreal = 1; 
options.issym = 1; 
[U,D] = eigs(B,d+2,0,options);  
lambda = diag(D);
[lambda_s,J] = sort(abs(lambda));
U = U(:,J); lambda = lambda(J);
T = U(:,2:d+1)';


% --- Creates and returns a handle to the GUI figure. 
function h1 = mani_LayoutFcn(policy)
% policy - create a new figure or use a singleton. 'new' or 'reuse'.

persistent hsingleton;
if strcmpi(policy, 'reuse') & ishandle(hsingleton)
    h1 = hsingleton;
    return;
end