1 | function [Y,details]=sde(DD,K,varargin) |
---|
2 | % [Y,details]=sde(DD,K,pars) |
---|
3 | % |
---|
4 | % |
---|
5 | % DD SQUARED distance matrix of the input vectors (e.g. euclidean distances) |
---|
6 | % |
---|
7 | % Optional: |
---|
8 | % |
---|
9 | % K number of neighbors |
---|
10 | % |
---|
11 | % pars Parameters |
---|
12 | % |
---|
13 | % pars.solver chooses the SDE solver: |
---|
14 | % pars.solver=0 CSDP (default) |
---|
15 | % pars.solver=1 SeDuMi |
---|
16 | % pars.solver=2 SDPT3 (very fast - not fully tested) |
---|
17 | % |
---|
18 | % pars.slack switches slack variables on or off |
---|
19 | % pars.slack=0 (default) |
---|
20 | % pars.slack=1 allows distances only to grow |
---|
21 | % pars.slack=2 allows both growing and shrinking of distances |
---|
22 | % pars.slack=3 allows distances only to shrink |
---|
23 | % (in combination with pars.factor=0 leads to |
---|
24 | % inequality) |
---|
25 | % |
---|
26 | % pars.factor Weight of penalty of Slack variables (only if pars.slac>0) |
---|
27 | % if pars.factor=0 the equality constraints become |
---|
28 | % inequalities (direction dependend on pars.slack) |
---|
29 | % default: pars.factor=0.9999 |
---|
30 | % |
---|
31 | % pars.inequ=0 default |
---|
32 | % pars.inequ=1 equivalent to pars.slack=3; pars.factor=0; |
---|
33 | % allos all distances to shrink by setting distance |
---|
34 | % preserving constraints to inequalities |
---|
35 | % |
---|
36 | % pars.verify Automatically increases K if graph is not connected |
---|
37 | % (Note: If no K is specified pars.verify=1) |
---|
38 | % pars.verify=1 (default) on |
---|
39 | % pars.verify=0 off |
---|
40 | % |
---|
41 | % pars.angles=1 preserves angles in local neighborhood (*DEFAULT*) |
---|
42 | % pars.angles=0 does not preserve angles in local neighborhood |
---|
43 | % |
---|
44 | % pars.repell Mx2 list of two points that should repell each other |
---|
45 | % pars.repell=[] (*default*) |
---|
46 | % |
---|
47 | % pars.rweight factor between 0 and 1, of how much the repelling |
---|
48 | % should weight compared to the normal objective |
---|
49 | % function |
---|
50 | % pars.rweight=0.3 (*default*) |
---|
51 | % |
---|
52 | % pars.kernel (independent of value) adds the original kernel matrix to the output |
---|
53 | % (details.K) |
---|
54 | % |
---|
55 | % Output: |
---|
56 | % |
---|
57 | % Y Resulting low dimensional vectors |
---|
58 | % (Rows are dimensions i.e. Y(1:d,:) is d-dimensional embedding |
---|
59 | % |
---|
60 | % details Includes all important details |
---|
61 | % k size of neighborhood |
---|
62 | % D all eigenvalues (ascending) |
---|
63 | % info SDP solver output |
---|
64 | % pars parameters |
---|
65 | % dual solution to dual problem |
---|
66 | % inder order of distance constraints (only interesting for |
---|
67 | % analysis of slack variables) |
---|
68 | % slack value of all slack variables |
---|
69 | % |
---|
70 | % |
---|
71 | % (TIP: From version 1.3 on, it is also possible to pass parameters |
---|
72 | % directly as strings |
---|
73 | % e.g. [Y,Det]=sde(D,3,'maxiter',20,'solver',0,'inequ'); |
---|
74 | % is equivalent to |
---|
75 | % pars.maxiter=20 |
---|
76 | % pars.solver=0 |
---|
77 | % pars.inequ=1; % (if the value is not specified, it is set to 1) |
---|
78 | % [Y,Det]=sde(D,3,pars); |
---|
79 | % |
---|
80 | % NOTE: sde requires a functionally copy of SeDuMi or CSDP in the path |
---|
81 | % (Download CSDP from http://www.nmt.edu/~borchers/csdp.html ) |
---|
82 | % (and SeDuMi from http://fewcal.kub.nl/sturm/software/sedumi.html ) |
---|
83 | % |
---|
84 | % Example: |
---|
85 | % tt=(linspace(0,1,100).^0.65).*3*pi+pi; |
---|
86 | % X=[tt.*cos(tt); tt.*sin(tt)]; % generates spiral input data |
---|
87 | % figure; |
---|
88 | % h1=scatter(X(1,:),X(2,:),140,tt,'filled'); % plots input data |
---|
89 | % title('INPUT DATA'); |
---|
90 | % set(h1,'MarkerEdgeColor',[0.5 0.5 0.5]); % draw edges around dots |
---|
91 | % drawnow; |
---|
92 | % axis equal; |
---|
93 | % [Y Det]=sde(distance(X),3,'maxiter',20,'inequ'); % runs SDE (with inequalities) |
---|
94 | % figure; |
---|
95 | % h2=scatter(Y(1,:),Y(2,:),140,tt,'filled'); % plots output data |
---|
96 | % set(h2,'MarkerEdgeColor',[0.5 0.5 0.5]); % draw edges around dots |
---|
97 | % title('OUTPUT DATA'); |
---|
98 | % axis equal; |
---|
99 | % |
---|
100 | % |
---|
101 | % You can execute this demo with |
---|
102 | % sde('demo'); |
---|
103 | % |
---|
104 | % (Version 1.3) |
---|
105 | % (copyright 2004 by Kilian Q. Weinberger: |
---|
106 | % http://www.seas.upenn.edu/~kilianw ) |
---|
107 | |
---|
108 | if(nargin==0) |
---|
109 | help sde; |
---|
110 | s=input('\n\nShall I run the SDE demo?','s'); |
---|
111 | if(length(s)>0 & isequal(upper(s(1)),'Y')) sde('demo');end; |
---|
112 | return; |
---|
113 | end; |
---|
114 | |
---|
115 | if(isequal(DD,'demo')) |
---|
116 | tt=(linspace(0,1,100).^0.65).*3*pi+pi; |
---|
117 | X=[tt.*cos(tt); tt.*sin(tt)]; % generates spiral input data |
---|
118 | figure; |
---|
119 | h1=scatter(X(1,:),X(2,:),140,tt,'filled'); % plots input data |
---|
120 | title('INPUT DATA'); |
---|
121 | set(h1,'MarkerEdgeColor',[0.5 0.5 0.5]); % draw edges around dots |
---|
122 | drawnow; |
---|
123 | axis equal; |
---|
124 | [Y Det]=sde(distance(X),3,'maxiter',20,'inequ'); % runs SDE (with inequalities) |
---|
125 | figure; |
---|
126 | h2=scatter(Y(1,:),Y(2,:),140,tt,'filled'); % plots output data |
---|
127 | set(h2,'MarkerEdgeColor',[0.5 0.5 0.5]); % draw edges around dots |
---|
128 | title('OUTPUT DATA'); |
---|
129 | axis equal; |
---|
130 | return; |
---|
131 | end; |
---|
132 | |
---|
133 | |
---|
134 | |
---|
135 | % extract variables |
---|
136 | if(length(varargin)>=1) |
---|
137 | if(~isstr(varargin{1})) |
---|
138 | pars=varargin{1}; |
---|
139 | for j=1:length(varargin)-1 |
---|
140 | varargin{j}=varargin{j+1}; |
---|
141 | end; |
---|
142 | end; |
---|
143 | |
---|
144 | for i=1:nargin-2 |
---|
145 | if(isstr(varargin{i})) |
---|
146 | if(i+1>nargin-2 | isstr(varargin{i+1})) val=1;else val=varargin{i+1};end; |
---|
147 | eval(['pars.' varargin{i} '=' sprintf('%i',val) ';']); |
---|
148 | end; |
---|
149 | end; |
---|
150 | end; |
---|
151 | |
---|
152 | |
---|
153 | [t1 t2]=size(DD); |
---|
154 | if(t1~=t2) |
---|
155 | % little check if distance matrix is really square |
---|
156 | error('First argument must be square(!!) distant matrix'); |
---|
157 | end; |
---|
158 | |
---|
159 | % Set standard parameters |
---|
160 | if(~exist('pars') | ~isfield(pars,'factor')) pars.factor=0.9999;end; |
---|
161 | if(~isfield(pars,'slack')) pars.slack=0;end; |
---|
162 | if(~isfield(pars,'verify')) pars.verify=1;end; |
---|
163 | if(~isfield(pars,'solver')) pars.solver=0;end; |
---|
164 | if(~isfield(pars,'angles')) pars.angles=1;end; |
---|
165 | if(~isfield(pars,'repell')) pars.repell=[];end; |
---|
166 | if(~isfield(pars,'rweight')) pars.rweight=0.3;end; |
---|
167 | if(~isfield(pars,'init')) pars.init=0;end; |
---|
168 | if(~isfield(pars,'maxiter')) pars.maxiter=50; end; |
---|
169 | if(~isfield(pars,'maxiter')) pars.maxiter=50; end; |
---|
170 | if(~isfield(pars,'inequ')) pars.inequ=0; end; |
---|
171 | |
---|
172 | |
---|
173 | |
---|
174 | % Calculate neighborhoods |
---|
175 | if(nargin<2) |
---|
176 | K=neighborsUDD(DD,3); |
---|
177 | else |
---|
178 | if(pars.verify) K=neighborsUDD(DD,K);end; |
---|
179 | end; |
---|
180 | |
---|
181 | if(pars.inequ==1) |
---|
182 | pars.slack=3; |
---|
183 | pars.factor=0; |
---|
184 | end; |
---|
185 | |
---|
186 | |
---|
187 | N=length(DD); |
---|
188 | [sorted,index] = sort(DD); |
---|
189 | neighbors = index(2:(1+K),:); |
---|
190 | |
---|
191 | if(isfield(pars,'ne')) |
---|
192 | neighbors=pars.ne; % not documented feature to pass as specific |
---|
193 | % neighborhood matrix as parameter |
---|
194 | end; |
---|
195 | |
---|
196 | |
---|
197 | % Constructing A and b - the constraints for the SDP |
---|
198 | nck=nchoosek(K+1,2); |
---|
199 | depth=N*nck*4; |
---|
200 | |
---|
201 | % Initialize matrix A |
---|
202 | A=zeros(depth,3); |
---|
203 | A(1:end,1)=vec(repmat(1:N*nck,4,1)); |
---|
204 | b=zeros(N*nck,1); |
---|
205 | |
---|
206 | % Initialize Matrix to check for redundant constraints |
---|
207 | DONE=sparse([],[],[],N,N,0); |
---|
208 | VALID=ones(N*nck,1); |
---|
209 | |
---|
210 | % just for detailed output, a mapping of input-points > constraints |
---|
211 | inder=zeros(depth,2); |
---|
212 | |
---|
213 | % KK = sdpvar(N,N); |
---|
214 | % F = set(KK>0) + set(sum(sum(KK))==0); |
---|
215 | % for i=1:N |
---|
216 | % ne = neighbors(:,i); |
---|
217 | % pairs=nchoosek([ne;i],2); |
---|
218 | % js=pairs(:,1); |
---|
219 | % ks=pairs(:,2); |
---|
220 | % for ii = 1:length(js) |
---|
221 | % F = F + set(KK(js(ii),js(ii))-2*KK(js(ii),ks(ii))+KK(ks(ii),ks(ii)) == DD(js(ii),ks(ii))); |
---|
222 | % end |
---|
223 | % end |
---|
224 | % solvesdp(F,-trace(KK)); |
---|
225 | % KK = double(KK) |
---|
226 | |
---|
227 | % Compute indicies |
---|
228 | |
---|
229 | |
---|
230 | if 1 |
---|
231 | allpairs = []; |
---|
232 | for i=1:N |
---|
233 | ne = neighbors(:,i); |
---|
234 | pairs=nchoosek([ne;i],2); |
---|
235 | allpairs = [allpairs;pairs]; |
---|
236 | end; |
---|
237 | allpairs = unique(sort(allpairs,2),'rows'); |
---|
238 | v1=sub2ind([N N],allpairs(:,2),allpairs(:,2)); |
---|
239 | v3=sub2ind([N N],allpairs(:,2),allpairs(:,1)); |
---|
240 | v4=sub2ind([N N],allpairs(:,1),allpairs(:,1)); |
---|
241 | % Setup problem |
---|
242 | X = sdpvar(N,N); |
---|
243 | s = sdpvar(length(allpairs),1); |
---|
244 | if(pars.slack==1) |
---|
245 | F = set(X >= 0) + set(s >= 0); |
---|
246 | F = F + set(sum(sum(X)) == 0) + set(X(v4) - 2*X(v3) + X(v1) + s == DD(v3)); |
---|
247 | obj = sum(s)*pars.factor-(1-pars.factor)*trace(X); |
---|
248 | else |
---|
249 | F = set(X >= 0); |
---|
250 | F = F + set(sum(sum(X)) == 0) + set(X(v4) - 2*X(v3) + X(v1) == DD(v3)); |
---|
251 | obj = -(1-pars.factor)*trace(X); |
---|
252 | end |
---|
253 | |
---|
254 | % Interpret in Primal cone |
---|
255 | [F,obj] = dualize(F,obj); |
---|
256 | |
---|
257 | % Solve |
---|
258 | %solvesdp(F,-obj,sdpsettings('solver','csdp','csdp.maxiter',50)); |
---|
259 | %X = double(X); |
---|
260 | end |
---|
261 | % |
---|
262 | Y=[]; |
---|
263 | details = []; |
---|
264 | return |
---|
265 | |
---|
266 | for i=1:N |
---|
267 | ne = neighbors(:,i); |
---|
268 | pairs=nchoosek([ne;i],2); |
---|
269 | |
---|
270 | js=pairs(:,1); |
---|
271 | ks=pairs(:,2); |
---|
272 | |
---|
273 | %find positions in output gram matrix |
---|
274 | v1=sub2ind([N N],ks,ks); |
---|
275 | v2=sub2ind([N N],js,ks); |
---|
276 | v3=sub2ind([N N],ks,js); |
---|
277 | v4=sub2ind([N N],js,js); |
---|
278 | |
---|
279 | |
---|
280 | pos=(i-1)*4*nck+1; |
---|
281 | |
---|
282 | vv=vec([v1 v2 v3 v4]'); |
---|
283 | % Eliminate doubles |
---|
284 | mask=DONE(v2)==0; |
---|
285 | |
---|
286 | if(pars.angles==0) |
---|
287 | mask=mask.*((js==i)+(ks==i)); |
---|
288 | end; |
---|
289 | |
---|
290 | VALID((i-1)*nck+1:(i-1)*nck+nck)=mask; |
---|
291 | DONE([v2;v3])=DONE([v2;v3])+1; |
---|
292 | A(pos:pos+4*nck-1,2)=vv; |
---|
293 | A(pos :4:pos+nck*4-4,3)=A(pos :4:pos+nck*4-4,3)+1; |
---|
294 | A(pos+1:4:pos+nck*4-3,3)=A(pos+1:4:pos+nck*4-3,3)-1; |
---|
295 | A(pos+2:4:pos+nck*4-2,3)=A(pos+2:4:pos+nck*4-2,3)-1; |
---|
296 | A(pos+3:4:pos+nck*4-1,3)=A(pos+3:4:pos+nck*4-1,3)+1; |
---|
297 | |
---|
298 | |
---|
299 | b(((i-1)*nck+1):(i*nck))= DD(sub2ind([N N],js,ks))'; |
---|
300 | |
---|
301 | inder((i-1)*nck+1:(i*nck),:)=[js ks]; |
---|
302 | end; |
---|
303 | A=sparse(A(:,1),A(:,2),A(:,3)); |
---|
304 | i=find(VALID>0); |
---|
305 | A=A(i,:); |
---|
306 | b=b(i,:); |
---|
307 | inder=inder(i,:); |
---|
308 | fprintf('Redundancy:%d%%\n',round(100-full(length(i))/(length(VALID))*100)); |
---|
309 | |
---|
310 | |
---|
311 | %% Add constraint to center points |
---|
312 | A=[ones(1,N^2);A]; |
---|
313 | b=[0;b]; |
---|
314 | i=[1;i]; |
---|
315 | |
---|
316 | |
---|
317 | |
---|
318 | if(isfield(pars,'c')) |
---|
319 | fprintf('Hidden Feature: Objective Function forced\n'); |
---|
320 | c=pars.c;%not documented feature to use weighted |
---|
321 | % objective function |
---|
322 | else |
---|
323 | c=0-vec(eye(N)); |
---|
324 | end; |
---|
325 | |
---|
326 | fprintf('Size of A: %ix%i \n\n\n',size(A)); |
---|
327 | |
---|
328 | flags.s=N; |
---|
329 | flags.l=0; |
---|
330 | |
---|
331 | %Y=[]; |
---|
332 | %details = []; |
---|
333 | %return |
---|
334 | |
---|
335 | %% Add repelling forces |
---|
336 | if(length(pars.repell>0)) |
---|
337 | [rp,temp]=size(pars.repell); |
---|
338 | c2=zeros(length(c),1); |
---|
339 | js=pars.repell(:,1); |
---|
340 | ks=pars.repell(:,2); |
---|
341 | v1=sub2ind([N N],ks,ks); |
---|
342 | v2=sub2ind([N N],js,ks); |
---|
343 | v3=sub2ind([N N],ks,js); |
---|
344 | v4=sub2ind([N N],js,js); |
---|
345 | |
---|
346 | for counter=1:rp |
---|
347 | c2(v1(counter))=1; |
---|
348 | c2(v2(counter))=-1;% |
---|
349 | c2(v3(counter))=-1;% |
---|
350 | c2(v4(counter))=1; |
---|
351 | end; |
---|
352 | |
---|
353 | c=c.*(1-pars.rweight)-c2.*pars.rweight; |
---|
354 | end; |
---|
355 | |
---|
356 | |
---|
357 | if(pars.slack==3) |
---|
358 | % Add Additive Slack Variables |
---|
359 | flags.l=length(i)-1; |
---|
360 | c=[ones(flags.l,1).*pars.factor;c.*(1-pars.factor)]; |
---|
361 | A=[[zeros(1,flags.l); speye(flags.l)] A]; |
---|
362 | end; |
---|
363 | |
---|
364 | |
---|
365 | if(pars.slack==2) |
---|
366 | % Add Additive Slack Variables (both additive and sub |
---|
367 | flags.l=(length(i)-1)*2; |
---|
368 | c=[ones(flags.l,1).*pars.factor;c.*(1-pars.factor)]; |
---|
369 | A=[[zeros(1,length(i)-1); speye(length(i)-1)] [zeros(1,length(i)-1); -speye(length(i)-1)] A]; |
---|
370 | end; |
---|
371 | |
---|
372 | if(pars.slack==1) |
---|
373 | % Add Additive Slack Variables |
---|
374 | flags.l=length(i)-1; |
---|
375 | c=[ones(flags.l,1).*pars.factor;c.*(1-pars.factor)]; |
---|
376 | A=[[zeros(1,flags.l); -speye(flags.l)] A]; |
---|
377 | end; |
---|
378 | |
---|
379 | if(pars.slack & 1==0) |
---|
380 | % Add Multiplicative Slack Variables |
---|
381 | flags.l=length(i)-1; |
---|
382 | c=[ones(flags.l,1).*pars.factor;c.*(1-pars.factor)]; |
---|
383 | A=[[zeros(1,flags.l); -sparse(spdiags(b(2:end),0,length(b)-1,length(b)-1))] A]; |
---|
384 | end; |
---|
385 | |
---|
386 | |
---|
387 | Kneigh=K; |
---|
388 | |
---|
389 | % Y=[]; |
---|
390 | % details=[]; |
---|
391 | % return |
---|
392 | |
---|
393 | switch(pars.solver) |
---|
394 | case{0}, |
---|
395 | OPTIONS.maxiter=50;%pars.maxiter; |
---|
396 | [x d z info]=csdp(A,b,c,flags,OPTIONS); |
---|
397 | K=mat(x(flags.l+1:flags.l+flags.s^2)); |
---|
398 | slack=x(1:flags.l); |
---|
399 | case{1}, |
---|
400 | OPTIONS.maxiter=pars.maxiter; |
---|
401 | [k d info]=sedumi(A,b,c,flags,OPTIONS); |
---|
402 | K=mat(k(end-flags.s^2+1:end)); |
---|
403 | slack=k(1:end-flags.s^2); |
---|
404 | %keyboard; |
---|
405 | |
---|
406 | case{2}, |
---|
407 | fprintf('converting into SQLP...'); |
---|
408 | [blk,A,c,b]=convert_sedumi(A',b,c,flags); |
---|
409 | fprintf('DONE\n'); |
---|
410 | OPTIONS.vers=1; |
---|
411 | OPTIONS.gam=0; |
---|
412 | OPTIONS.predcorr=1; |
---|
413 | OPTIONS.expon=[1 1]; |
---|
414 | OPTIONS.gaptol=1.0000e-08; |
---|
415 | OPTIONS.inftol=1.0000e-08; |
---|
416 | OPTIONS.steptol=1.0000e-06; |
---|
417 | OPTIONS.maxit=pars.maxiter; |
---|
418 | OPTIONS.printyes=1; |
---|
419 | OPTIONS.scale_data=0; |
---|
420 | OPTIONS.randnstate=0; |
---|
421 | OPTIONS.spdensity=0.5000; |
---|
422 | OPTIONS.rmdepconstr=0; |
---|
423 | OPTIONS.cachesize=256; |
---|
424 | OPTIONS.smallblkdim=15; |
---|
425 | |
---|
426 | if(pars.init==1) |
---|
427 | fprintf('Initializing ...'); |
---|
428 | % Initialize with input points |
---|
429 | [X0,y0,Z0]=infeaspt(blk,A',c,b); |
---|
430 | X0{length(X0)}=getgram(DD)+eye(N).*0.001; |
---|
431 | fprintf('DONE\n'); |
---|
432 | [obj,Kt,d,z]=sqlp(blk,A,c,b,{},X0,y0,Z0); |
---|
433 | else |
---|
434 | OPTIONS.maxit=pars.maxiter; |
---|
435 | [obj,Kt,d,z]=sqlp(blk,A,c,b,OPTIONS); |
---|
436 | end; |
---|
437 | K=Kt{length(Kt)}; |
---|
438 | info=0; |
---|
439 | if(pars.slack)slack=Kt{1}(:);end; |
---|
440 | end; |
---|
441 | %keyboard; |
---|
442 | [V D]=eig(K); |
---|
443 | |
---|
444 | V=V*sqrt(D); |
---|
445 | D=diag(D); |
---|
446 | |
---|
447 | Y=(V(:,end:-1:1))'; |
---|
448 | |
---|
449 | if(isfield(pars,'kernel')) |
---|
450 | details.K=K; |
---|
451 | end; |
---|
452 | details.k=Kneigh; |
---|
453 | details.D=D; |
---|
454 | details.info=info; |
---|
455 | details.pars=pars; |
---|
456 | details.dual=d; |
---|
457 | details.inder=inder; |
---|
458 | |
---|
459 | |
---|
460 | if(pars.slack)details.slack=slack;end; |
---|
461 | |
---|
462 | return; |
---|
463 | |
---|
464 | |
---|
465 | |
---|
466 | |
---|
467 | |
---|
468 | |
---|
469 | |
---|
470 | |
---|
471 | |
---|
472 | |
---|
473 | |
---|
474 | |
---|
475 | |
---|
476 | |
---|
477 | |
---|
478 | |
---|
479 | |
---|
480 | |
---|
481 | |
---|
482 | %%% SIMPLE BUT USEFUL UTIL FUNCTIONS |
---|
483 | |
---|
484 | function neighbors=getneighborsUDD(DD,K); |
---|
485 | ne=getneighborsDD(DD,K); |
---|
486 | N=length(DD); |
---|
487 | |
---|
488 | for i=1:N |
---|
489 | neighbors{i}=[]; |
---|
490 | end; |
---|
491 | |
---|
492 | for i=1:N |
---|
493 | neighbors{i}=merge(sort(neighbors{i}),sort(ne(:,i))); |
---|
494 | for j=1:K |
---|
495 | neighbors{ne(j,i)}=merge(neighbors{ne(j,i)}, i); |
---|
496 | end; |
---|
497 | end; |
---|
498 | |
---|
499 | |
---|
500 | |
---|
501 | function result=merge(x,y); |
---|
502 | result=unique([x(:);y(:)]); |
---|
503 | |
---|
504 | |
---|
505 | function v=vec(M); |
---|
506 | v=M(:); |
---|
507 | |
---|
508 | |
---|
509 | function k=neighborsUDD(DD,K); |
---|
510 | N=length(DD); |
---|
511 | if(nargin<2) |
---|
512 | K=2; |
---|
513 | end; |
---|
514 | k=K; |
---|
515 | while(k<N & (1-connectedUDD(DD,k))) |
---|
516 | k=k+1; |
---|
517 | fprintf('Trying K=%i \n',k); |
---|
518 | end; |
---|
519 | |
---|
520 | |
---|
521 | |
---|
522 | function result=connectedUDD(DD,K); |
---|
523 | % result = connecteddfs (X,K) |
---|
524 | % |
---|
525 | % X input vector |
---|
526 | % K number of neighbors |
---|
527 | % |
---|
528 | % Returns: result = 1 if connected 0 if not connected |
---|
529 | |
---|
530 | if(nargin<2) |
---|
531 | fprintf('Number of Neighbors not specified!\nSetting K=4\n'); |
---|
532 | K=4; |
---|
533 | end; |
---|
534 | N=length(DD); |
---|
535 | |
---|
536 | |
---|
537 | ne=getneighborsUDD(DD,K); |
---|
538 | |
---|
539 | |
---|
540 | maxSize=0; |
---|
541 | for i=1:N |
---|
542 | if(length(ne{i})>maxSize) maxSize=length(ne{i});end; |
---|
543 | end; |
---|
544 | neighbors=ones(maxSize,N); |
---|
545 | for i=1:N |
---|
546 | neighbors(1:length(ne{i}),i)=ne{i}; |
---|
547 | end; |
---|
548 | oldchecked=[]; |
---|
549 | checked=[1]; |
---|
550 | |
---|
551 | while((size(checked)<N) & (length(oldchecked)~=length(checked))) |
---|
552 | next=neighbors(:,checked); |
---|
553 | next=transpose(sort(next(:))); |
---|
554 | next2=[next(2:end) 0]; |
---|
555 | k=find(next-next2~=0); |
---|
556 | next=next(k); |
---|
557 | |
---|
558 | oldchecked=checked; |
---|
559 | checked=neighbors(:,next(:)); |
---|
560 | checked=transpose(sort(checked(:))); |
---|
561 | checked2=[checked(2:end) 0]; |
---|
562 | k=find(checked-checked2~=0); |
---|
563 | checked=checked(k); |
---|
564 | end; |
---|
565 | result=(length(checked)==N); |
---|
566 | |
---|
567 | |
---|
568 | |
---|
569 | function X = mat(x,n) |
---|
570 | if nargin < 2 |
---|
571 | n = floor(sqrt(length(x))); |
---|
572 | if (n*n) ~= length(x) |
---|
573 | error('Argument X has to be a square matrix') |
---|
574 | end |
---|
575 | end |
---|
576 | X = reshape(x,n,n); |
---|
577 | |
---|
578 | |
---|
579 | |
---|
580 | |
---|
581 | function neighbors=getneighborsDD(DD,K); |
---|
582 | N=length(DD); |
---|
583 | [sorted,index] = sort(DD); |
---|
584 | neighbors = index(2:(1+K),:); |
---|