[37] | 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),:); |
---|