function [ newX,A ] = CDA( y,n_fcn,n_size,alpha,echo,error )
%CDA Summary of this function goes here
% Detailed explanation goes here
D=L2_distance(y',y');
N = size(D,1);
INF = 1000*max(max(D))*N;
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');
end
elseif n_fcn=='epsilon'
epsilon = n_size;
else
error('Neighborhood function must be either epsilon or k');
end
if n_fcn == 'k'
[tmp, ind] = sort(D);
for i=1:N
D(i,ind((2+K):end,i)) = INF;
end
elseif n_fcn == 'epsilon'
warning off %% Next line causes an unnecessary warning, so turn it off
D = D./(D<=epsilon);
D = min(D,INF);
warning on
end
D = min(D,D'); %% Make sure distance matrix is symmetric
% if (overlay == 1)
E = int8(1-(D==INF)); %% Edge information for subsequent graph overlay
% end
% 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']);
% end
end
Y_GD=D;
%输入数据y的图距离
[my,ny]=size(y);
N=my;
train_len=echo*my;
P=round(10*rand(3,N))';
%初始化降维目标数据
Ecca=inf;
k=1;
% lambda_temp=lambda;
alpha_temp=alpha;
P_ED=squareform(pdist(P));
pi=1;
% Ecca_miu=zeros(N,N);
% Ecca_p=zeros(N,N);
A=[];
for i=1:N
lampda(i)=max(P_ED(i,:));
end
while (k<=train_len&&error<Ecca)
E_temp=0;
Ecca_miu=zeros(N,N);
Ecca_p=zeros(N,N);
for i=1:N
% lampda(i)=max(P_ED(i,:));
for j=1:N
if i~=j
if P_ED(i,j)>Y_GD(i,j)
if lampda(i)>=P_ED(i,j)
delta_Ecca_miu=2*(Y_GD(i,j)-P_ED(i,j))*1*(P(i,:)-P(j,:))/P_ED(i,j);
else
delta_Ecca_miu=0;
end
P(j,:)= P(j,:)+alpha_temp(k)*delta_Ecca_miu;
Ecca_miu(i,j)=0.5*(Y_GD(i,j)-P_ED(i,j))^2*1;
% Ecca_miu(i,j)=Ecca_temp;
% P(j,:)= P(j,:)+alpha_temp(k)*delta_Ecca_miu;
else %P_ED(i,j)<=Y_GD(i,j)
if lampda(i)>=P_ED(i,j)
delta_Ecca_p=-4*(Y_GD(i,j)^2-P_ED(i,j)^2)*1*(P(i,:)-P(j,:));
else
delta_Ecca_p=0;
end
% Ecca_p(i,j)=Ecca_temp;
P(j,:)= P(j,:)-alpha_temp(k)*1/(4*Y_GD(i,j)^2)*delta_Ecca_p;
Ecca_p(i,j)=0.5*1/(4*Y_GD(i,j)^2)*(Y_GD(i,j)^2-P_ED(i,j)^2)^2*1;
end
else
Ecca_miu(i,j)=0;
Ecca_p(i,j)=0;
end
lampda(i)= lampda(i)*(pi/(n_size/N))^(1/3);
end
end
for i=1:N
for j=1:N
if i~=j
if P_ED(i,j)>Y_GD(i,j)
E_temp=E_temp+Ecca_miu(i,j);
else
E_temp=E_temp+1/(4*Y_GD(i,j)^2)*Ecca_p(i,j);
end
if P_ED(i,j)>Y_GD(i,j)
E_temp=E_temp+Ecca_miu(i,j);
else
E_temp=E_temp+1/(4*Y_GD(i,j)^2)*Ecca_p(i,j);
end
% else
% E_temp=0;
end
end
end
Ecca=E_temp;
k=k+1;
alpha_temp=potency_curve(alpha,alpha/100,k);
P_ED=squareform(pdist(P));
for i=1:N
lampda(i)=max(P_ED(i,:));
end
A(k)=Ecca;
end
newX=P;
% plot(A);
% axis([0,200]);
% axis auto;
function vals = potency_curve(v0,vn,l)
vals = v0 * (vn/v0).^([0:(l-1)]/(l-1));
评论0