ID:46063
Title:wire through 02
Author:Alfonso Nieto-Castanon
Date:2008-05-02 00:43:22
Score:15619.7608
Result:155461.00 (cyc: 30, node: 3262)
CPU Time:48.4023
Status:Passed
Comments:
Based on:none
Code:
function W=solver(B);

totEE=0;
W=zeros(0,4);
sB=size(B);
B0=zeros(sB);
idx0=find(B(:));
[b,I,J]=unique(B(idx0));

L0=zeros(size(b));
for n1=1:length(b),
	idx10{n1}=find(J==n1);
	L0(n1)=length(idx10{n1});
end
L0(L0==1)=0;
[nill,sortN]=sort(L0.*b);

for nsortN=length(sortN):-1:1+sum(L0==0),
	N=sortN(nsortN);
	L1=L0(N);
	idx1=idx10{N};

	if L1>1,
		R=cell(2*L1-1,2*L1-1);C=R;Q=R; % connecting paths
		D=zeros(2*L1-1,2*L1-1); D(:)=inf; % distances
		E=zeros(2*L1-1,2*L1-1); % gain
		IDX=zeros(2*L1-1,1); IDX(1:L1)=1; % valid clusters
		[r,c]=ind2sub(sB,idx0(idx1));
		for n1=1:L1, RR{n1}=r(n1); CC{n1}=c(n1); QQ{n1}=[]; end % clusters
		EE=zeros(2*L1-1,1); EE(1:L1)=b(N);
		WW=cell(2*L1-1,1);
		
		% first pass
		LL=2/L1/(L1-1);
		e=0; idx=[];
		for n1=1:L1, 
			for n2=n1+1:L1, 
				[R{n1,n2},C{n1,n2},Q{n1,n2},D(n1,n2)]=findpath(B,r(n1),c(n1),r(n2),c(n2)); 
				E(n1,n2)=EE(n1)+EE(n2)-D(n1,n2);
				if E(n1,n2)>e || (E(n1,n2)==e&&e>0&&rand<LL), e=E(n1,n2); idx=[n1,n2]; end
			end
		end
		K=L1; 
		if e>0, % something
			while ~isempty(idx),
				% new cluster
				K=K+1; 
				Rnow=R{idx(1),idx(2)}; Cnow=C{idx(1),idx(2)}; Qnow=Q{idx(1),idx(2)}; 
				RR{K}=cat(1,RR{idx(1)},Rnow((2:end-1)'),RR{idx(2)});
				CC{K}=cat(1,CC{idx(1)},Cnow((2:end-1)'),CC{idx(2)});
				QQ{K}=cat(1,QQ{idx(1)},Qnow,QQ{idx(2)});
				EE(K)=EE(idx(1))+EE(idx(2))-D(idx(1),idx(2));
				WW{K}=cat(1,WW{idx(1)},WW{idx(2)},[Rnow(1:end-1),Cnow(1:end-1),Rnow(2:end),Cnow(2:end)]);
				IDX(idx)=0;
				n1=K; N2=find(IDX);
				IDX(K)=1;
				if ~isempty(N2), 
					for n2=N2(:)', 
						if idx(1)<n2, mi1=idx(1);ma1=n2; else, mi1=n2;ma1=idx(1); end
						if idx(2)<n2, mi2=idx(2);ma2=n2; else, mi2=n2;ma2=idx(2); end
						if D(mi1,ma1)<D(mi2,ma2)
							R{n2,n1}=R{mi1,ma1};
							C{n2,n1}=C{mi1,ma1};
							Q{n2,n1}=Q{mi1,ma1};
							D(n2,n1)=D(mi1,ma1);
							E(n2,n1)=EE(n1)+EE(n2)-D(n2,n1);
						else, 
							R{n2,n1}=R{mi2,ma2};
							C{n2,n1}=C{mi2,ma2};
							Q{n2,n1}=Q{mi2,ma2};
							D(n2,n1)=D(mi2,ma2);
							E(n2,n1)=EE(n1)+EE(n2)-D(n2,n1);
						end
						if D(idx(1),idx(2))>1,
							% recompute distances
							LL=1/length(N2);
							[Rt,Ct,Qt,Dt]=findpath(B,Rnow((2:end-1)'),Cnow((2:end-1)'),RR{n2},CC{n2},1); 
							%disp([Rt,Ct]);
							if Dt<=D(n2,n1), R{n2,n1}=Rt;C{n2,n1}=Ct;Q{n2,n1}=Qt;D(n2,n1)=Dt; end
							E(n2,n1)=EE(n1)+EE(n2)-D(n2,n1);
						end
					end
					
					idx=[];
					IDX1=find(IDX);N1=length(IDX1);d=D(IDX1,IDX1);
					[sortd,idxd]=sort(d(:)); % closest path
					[r1,c1]=ind2sub([N1,N1],idxd(1:N1*(N1-1)/2));
					idxfirst=find(E(IDX1(r1)+size(E,1)*(IDX1(c1)-1))>max(EE(IDX1(r1)),EE(IDX1(c1))));  % first increasing gain
					if ~isempty(idxfirst), idx=[IDX1(r1(idxfirst(1))),IDX1(c1(idxfirst(1)))]; end
				else, idx=[]; end
			end
		end
		if K>L1,
			%for n1=1:K,EE(n1)=EE(n1)-25*sum(B0(RR{n1}+sB(1)*(CC{n1}-1))>0);end
			[maxEE,K]=max(EE(L1+1:end)); K=L1+K;
			if maxEE>0 && length(RR{K})>1,
				totEE=totEE+maxEE; 
				%disp(['On node ',num2str(b(N)),' score ',num2str(maxEE)]);
				W=cat(1,W,WW{K});
				for n3=1:length(RR{K}), 
					if ~B(RR{K}(n3),CC{K}(n3)), B(RR{K}(n3),CC{K}(n3))=-1; end; % paths
					if B0(RR{K}(n3),CC{K}(n3))>0, 
						W=cat(1,W,[RR{K}(n3),CC{K}(n3),RR{K}(n3),CC{K}(n3)]);
					end
				end
				for n3=1:size(QQ{K},1), 
					B(QQ{K}(n3,1),QQ{K}(n3,2))=1; % corners
				end
				for n3=1:length(RR{K}), 
					B0(RR{K}(n3),CC{K}(n3))=1; % paths
				end
			end
		end
	end
end
%disp(totEE);
%figure(99); imagesc([B,B0]);

function [R,C,Q,L]=findpath(B,r1,c1,r2,c2,xtra);
% finds straight path from r1 c1 to r2 c2 in map B
R=[];C=[];Q=[];L=inf;
L1=length(r1); L2=length(r2);
A=abs(r1(:,ones(1,L2))-r2(:,ones(1,L1))')+abs(c1(:,ones(1,L2))-c2(:,ones(1,L1))');
[a,idx]=sort(A(:));
for n1=1:length(a),
	if L<=a(n1),return;end
	i2=ceil(idx(n1)/L1);
	i1=idx(n1)-(i2-1)*L1;
	% From B(c1(i1),r1(i1)) to B(c2(i2),r2(i2))
	if c1(i1)==c2(i2), % vertical
		C0=c1(i1);
		if r2(i2)<r1(i1), R0=(r1(i1):-1:r2(i2))';
		else, R0=(r1(i1):r2(i2))'; end
		L0=length(R0)-1;
		if L0<2||~any(B(R0(2:end-1),C0)), % clean straight
			R=R0;
			C=C0(ones(L0+1,1));
			Q=[r1(i1),c1(i1);r2(i2),c2(i2)];
			L=L0;
			return;
		end
	end
	if r1(i1)==r2(i2), % horizontal
		R0=r1(i1);
		if c2(i2)<c1(i1), C0=(c1(i1):-1:c2(i2))';
		else, C0=(c1(i1):c2(i2))'; end
		L0=length(C0)-1;
		if L0<2||~any(B(R0,C0(2:end-1))), % clean straight
			C=C0;
			R=R0(ones(L0+1,1));
			Q=[r1(i1),c1(i1);r2(i2),c2(i2)];
			L=L0;
			return;
		end
	end
	if c1(i1)<c2(i2), minc=c1(i1);maxc=c2(i2); idxc=1; else, minc=c2(i2);maxc=c1(i1); idxc=2; end
	b1=B(:,c1(i1)); b1(r1(i1))=0; 
	bb1a=cumsum(b1>0);
	bb1b=bb1a; bb1b(1)=bb1a(end);bb1b(2:end)=bb1a(end)-bb1a(1:end-1); 
	bb1=max(0,bb1a-bb1a(r1(i1)))+max(0,bb1b-bb1b(r1(i1))); 
	d1a=cumsum(b1~=0);
	d1b=d1a; d1b(1)=d1a(end);d1b(2:end)=d1a(end)-d1a(1:end-1); 
	d1=max(0,d1a-d1a(r1(i1)))+max(0,d1b-d1b(r1(i1)));
	b2=B(:,c2(i2)); b2(r2(i2))=0; 
	bb2a=cumsum(b2>0);
	bb2b=bb2a; bb2b(1)=bb2a(end);bb2b(2:end)=bb2a(end)-bb2a(1:end-1); 
	bb2=max(0,bb2a-bb2a(r2(i2)))+max(0,bb2b-bb2b(r2(i2))); 
	d2a=cumsum(b2~=0);
	d2b=d2a; d2b(1)=d2a(end);d2b(2:end)=d2a(end)-d2a(1:end-1); 
	d2=max(0,d2a-d2a(r2(i2)))+max(0,d2b-d2b(r2(i2)));
	validcorner1=find(b1==0&b2==0&bb1==0&bb2==0);
	if ~isempty(validcorner1),
		e=any(B(validcorner1,minc+1:maxc-1)>0,2);
		validcorner2=find(e==0);
		if ~isempty(validcorner2);
			f=sum(B(validcorner1(validcorner2),minc+1:maxc-1)~=0,2);
			validcorner=validcorner1(validcorner2);
			d=abs(validcorner-r1(i1))+abs(c2(i2)-c1(i1))+abs(validcorner-r2(i2))+25*(d1(validcorner)+d2(validcorner)+f);
			[md,id]=min(d);
			if md<L,
				if idxc==1, C=cat(1,c1(i1)+zeros(abs(validcorner(id)-r1(i1)),1),(c1(i1):c2(i2))',c2(i2)+zeros(abs(validcorner(id)-r2(i2)),1));
				else, 		C=cat(1,c1(i1)+zeros(abs(validcorner(id)-r1(i1)),1),(c1(i1):-1:c2(i2))',c2(i2)+zeros(abs(validcorner(id)-r2(i2)),1)); end
				if validcorner(id)>=r1(i1),
					if r2(i2)>=validcorner(id),R=cat(1,(r1(i1):validcorner(id))',validcorner(id)+zeros(abs(c2(i2)-c1(i1)),1),(validcorner(id)+1:r2(i2))');
					else, 					   R=cat(1,(r1(i1):validcorner(id))',validcorner(id)+zeros(abs(c2(i2)-c1(i1)),1),(validcorner(id)-1:-1:r2(i2))'); end
				else,
					if r2(i2)>=validcorner(id),R=cat(1,(r1(i1):-1:validcorner(id))',validcorner(id)+zeros(abs(c2(i2)-c1(i1)),1),(validcorner(id)+1:r2(i2))');
					else, 					   R=cat(1,(r1(i1):-1:validcorner(id))',validcorner(id)+zeros(abs(c2(i2)-c1(i1)),1),(validcorner(id)-1:-1:r2(i2))'); end
				end
				Q=[r1(i1),c1(i1);r2(i2),c2(i2);validcorner(id),minc; validcorner(id),maxc];
				L=md;
				if L<=a(n1),return;end
			end
		end
	end
	if r1(i1)<r2(i2), minr=r1(i1);maxr=r2(i2); idxr=1; else, minr=r2(i2);maxr=r1(i1); idxr=2; end
	b1=B(r1(i1),:)'; b1(c1(i1))=0; 
	bb1a=cumsum(b1>0);
	bb1b=bb1a; bb1b(1)=bb1a(end);bb1b(2:end)=bb1a(end)-bb1a(1:end-1); 
	bb1=max(0,bb1a-bb1a(c1(i1)))+max(0,bb1b-bb1b(c1(i1))); 
	d1a=cumsum(b1~=0);
	d1b=d1a; d1b(1)=d1a(end);d1b(2:end)=d1a(end)-d1a(1:end-1); 
	d1=max(0,d1a-d1a(c1(i1)))+max(0,d1b-d1b(c1(i1)));
	b2=B(r2(i2),:)'; b2(c2(i2))=0; 
	bb2a=cumsum(b2>0);
	bb2b=bb2a; bb2b(1)=bb2a(end);bb2b(2:end)=bb2a(end)-bb2a(1:end-1); 
	bb2=max(0,bb2a-bb2a(c2(i2)))+max(0,bb2b-bb2b(c2(i2))); 
	d2a=cumsum(b2~=0);
	d2b=d2a; d2b(1)=d2a(end);d2b(2:end)=d2a(end)-d2a(1:end-1); 
	d2=max(0,d2a-d2a(c2(i2)))+max(0,d2b-d2b(c2(i2)));
	validcorner1=find(b1==0&b2==0&bb1==0&bb2==0);
	if ~isempty(validcorner1),
		e=any(B(minr+1:maxr-1,validcorner1)>0,1)';
		validcorner2=find(e==0);
		if ~isempty(validcorner2);
			f=sum(B(minr+1:maxr-1,validcorner1(validcorner2))~=0,1)';
			validcorner=validcorner1(validcorner2);
			d=abs(validcorner-c1(i1))+abs(r2(i2)-r1(i1))+abs(validcorner-c2(i2))+25*(d1(validcorner)+d2(validcorner)+f);
			[md,id]=min(d);
			if md<L,
				if idxr==1, R=cat(1,r1(i1)+zeros(abs(validcorner(id)-c1(i1)),1),(r1(i1):r2(i2))',r2(i2)+zeros(abs(validcorner(id)-c2(i2)),1));
				else, 		R=cat(1,r1(i1)+zeros(abs(validcorner(id)-c1(i1)),1),(r1(i1):-1:r2(i2))',r2(i2)+zeros(abs(validcorner(id)-c2(i2)),1)); end
				if validcorner(id)>=c1(i1),
					if c2(i2)>=validcorner(id),C=cat(1,(c1(i1):validcorner(id))',validcorner(id)+zeros(abs(r2(i2)-r1(i1)),1),(validcorner(id)+1:c2(i2))');
					else, 					   C=cat(1,(c1(i1):validcorner(id))',validcorner(id)+zeros(abs(r2(i2)-r1(i1)),1),(validcorner(id)-1:-1:c2(i2))'); end
				else,
					if c2(i2)>=validcorner(id),C=cat(1,(c1(i1):-1:validcorner(id))',validcorner(id)+zeros(abs(r2(i2)-r1(i1)),1),(validcorner(id)+1:c2(i2))');
					else, 					   C=cat(1,(c1(i1):-1:validcorner(id))',validcorner(id)+zeros(abs(r2(i2)-r1(i1)),1),(validcorner(id)-1:-1:c2(i2))'); end
				end
				Q=[r1(i1),c1(i1);r2(i2),c2(i2);minr,validcorner(id); maxr,validcorner(id)];
				L=md;
				if L<=a(n1),return;end
			end
		end
	end
end

if isinf(L),C=[];R=[];Q=[]; end
return