📄 mani.m
字号:
[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'));endfunction 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'));endfunction 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
h1 = figure(...
'Units','characters',...
'Color',[0.925490196078431 0.913725490196078 0.847058823529412],...
'Colormap',[0 0 0.5625;0 0 0.625;0 0 0.6875;0 0 0.75;0 0 0.8125;0 0 0.875;0 0 0.9375;0 0 1;0 0.0625 1;0 0.125 1;0 0.1875 1;0 0.25 1;0 0.3125 1;0 0.375 1;0 0.4375 1;0 0.5 1;0 0.5625 1;0 0.625 1;0 0.6875 1;0 0.75 1;0 0.8125 1;0 0.875 1;0 0.9375 1;0 1 1;0.0625 1 1;0.125 1 0.9375;0.1875 1 0.875;0.25 1 0.8125;0.3125 1 0.75;0.375 1 0.6875;0.4375 1 0.625;0.5 1 0.5625;0.5625 1 0.5;0.625 1 0.4375;0.6875 1 0.375;0.75 1 0.3125;0.8125 1 0.25;0.875 1 0.1875;0.9375 1 0.125;1 1 0.0625;1 1 0;1 0.9375 0;1 0.875 0;1 0.8125 0;1 0.75 0;1 0.6875 0;1 0.625 0;1 0.5625 0;1 0.5 0;1 0.4375 0;1 0.375 0;1 0.3125 0;1 0.25 0;1 0.1875 0;1 0.125 0;1 0.0625 0;1 0 0;0.9375 0 0;0.875 0 0;0.8125 0 0;0.75 0 0;0.6875 0 0;0.625 0 0;0.5625 0 0],...
'IntegerHandle','off',...
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -