function [u,val,X,newP,newg] = stemprocedure(G,dist,n,k,P,g,nmd,v,d)
%原图G，第n+1个点到其他的距离dist，n为现有的点数,k为k近邻的个数,P为先驱矩阵, g为测地线矩阵, nmd为特征值，v为特征向量，d为维数
neighbors = zeros(n+1,1);
for i = 1:n
G(n+1,i) = dist(i);
end
%得到k近邻
for i = 1:n+1
work = G(i,:);
[~,index] = sort(work,'ascend');
neighbors(i)= index(k+1);
end
%得到A D
A = zeros(n+1,2);
counter = 1;
for i = 1:n
delta1 = G(i,neighbors(i));
delta2 = G(n+1,i);
if(delta1>delta2)
A(counter,1) = i;
A(counter,2) = n+1;
counter = counter+1;
end
end
D = zeros(n+1,2);
counter = 1;
for i = 1:n
delta1 = G(i,neighbors(i));
delta2 = G(n+1,i);
if delta1>delta2
taoi = neighbors(i);
delta3 = G(taoi,i);
delta4 = G(neighbors(taoi),taoi);
if delta3 >delta4
D(counter,1) = i;
D(counter,2) = taoi;
counter = counter+1;
end
end
end
for i = 1:n
G(i,n+1) = dist(i);
end
%得到影响的点对
F = zeros(2,10000);
F(:,1) = [-1;-1];
counter = 1;% for F
for i = 1:size(D,2)
if D(i,1) ~= 0
[Fab,Fsize] = ConstructFab([D(i,1),D(i,2)],G(1:n,1:n),P);
for j = 1:Fsize
x = Fab(j,1);
y = Fab(j,2);
if ismember(x,F(1,:))
[~,loc] = ismember(x,F(1,:));
for k = 1:loc
if F(1,k) == x;
if y~=F(2,k)
counter = counter+1;
F(:,counter) = [x;y];
end
end
end
elseif ismember(x,F(2,:))
[~,loc] = ismember(x,F(2,:));
for k = 1:loc
if F(2,k) == x
if y~=F(1,k)
counter = counter+1;
F(:,counter) = [x;y];
end
end
end
else
counter = counter+1;
F(:,counter) = [x;y];
end
end
end
end
F = F(:,2:counter);
for i = 1:size(F,2)
s = F(1,i);
t = F(2,i);
G(s,t) = inf;
G(t,s) = inf;
end
%更新g矩阵
u = unique(F(1,:));
for i = 1:length(u)
[~,loc] = ismember(u(i),F(1,:));
Cu = zeros(1,10000);
counter = 1; %for Cu
for k = 1:loc
if F(1,k) == u(i)
Cu(counter) = F(2,k);
counter = counter+1;
end
end
g = modifieddijkstra(u(i),Cu(1:counter-1),g,G(1:n,1:n));
end
%记录Vn+1的前驱矩阵,g增维，P增维
[gn1,Tn1] = dijkstra(G(1:n+1,1:n+1),n+1);
newg = [g,zeros(n,1);gn1];
newP = [P,zeros(n,1);Tn1];
%填写前n个点到n+1个点的前驱
for i = 1:n
fore = newP(n+1,i);
if fore == n+1
newP(i,n+1) = i;
elseif fore~= 0
newP(i,n+1) = P(i,fore);
else
newP(i,n+1) = i;
end
end
%加入新边以后对于最短路径的更新
for i=1:n+1
if G(n+1,i)~=inf
for j=i:n+1
if G(n+1,j)~=inf && i~=j
if G(n+1,i) +G(n+1,j) < newg(i,j)
newg = UpdateInsert(i , j ,G,Tn1,newg,n+1);
end
end
end
end
end
%保持newg的对称性
newg(1:n,n+1) = newg(n+1,1:n);
f = zeros(1,n);
sumlj = 0;
%计算f
for j = 1:n
for l = 1:n
sumlj = sumlj + newg(j,l).*newg(j,l);
end
end
sumln1 = 0;
for l = 1:n
sumln1 = sumln1 + newg(l,n+1).*newg(l,n+1);
end
for i = 1:n
sumij = 0;
for j = 1:n
sumij = sumij + newg(i,j).*newg(i,j);
end
f(i) = 0.5*(sumij/n-sumlj/(n*n)+sumln1/n-newg(i,n+1).*newg(i,n+1));
% f(i) = 0.5*(sumij/n-newg(i,n+1).*newg(i,n+1));
end
%计算X,xn+1
X = zeros(d,n+1);
xn1 = zeros(1,d);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:d
X(i,1:n) = sqrt(nmd(i)).*real(v(:,i));
xn1(i) = (1/sqrt(nmd(i))).*real(v(:,i)')*f';%%%%%%%%出现问题的地方%%%%%%%%%%%%%%%%%%
%%%归一化
% xn1(i) = xn1(i)./10^5;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X(:,n+1) = xn1;
H = -1/(n+1)*ones(n+1);
H = H+diag(ones(1,n+1));
Bnew = -H*(newg.*newg)*H/2;
Z = real(Bnew)*X';
[~,V] = qr(Z);
Zxin = V'*real(Bnew)*V;
[u,val] = eig(real(Zxin));
[~,sorth] = sort(diag(val),'descend');
val = diag(val(sorth,sorth)); %返回特征值和特征向量
u = V*u;
u = u(:,sorth);
end
增量ISOMAP算法实现-源码-MATLAB
需积分: 34 29 浏览量
2015-04-01
14:57:56
上传
评论 2
收藏 8.62MB ZIP 举报
ISOMAP.zip （14个子文件） 
dijk.m 3KB batprocess.m 1KB dfun.m 85B stemprocedure.m 5KB swiss_roll_data.mat 782KB D.mat 7.94MB gudingjidian_bodongjiance.m 616B L2_distance.m 2KB modifieddijkstra.m 1KB UpdateInsert.m 2KB ConstructFab.m 1KB dijkstra.m 2KB OptimalOrder.m 1KB InitializeEdgeWeightIncrease.m 1KB资源评论
