Winner Paulo Uribe (turbf1)

2002-05-23 00:00:00 UTC

# n1wd1wt

by wc

Status: Passed
Results: Average strain = 2375.817
CPU Time: 114.585
Score: 8443.1
Submitted at: 2002-05-18 08:55:11 UTC
Scored at: 2002-05-18 09:01:12 UTC

Current Rank: 714th
Based on: n1wd1 (diff)
Basis for: n1wd1wt2 (diff)

Code
```function xy = solver(kd,bx)

npts = length(kd);
if npts == 1,
xy=bx([1 3]);
return
elseif npts == 2,
p0=bx([1 3]);
xy=p0;
p1=bx([2 4]);
l=kd(1,2);
if l == -1,
xy(2,:)=p1;
return
end
h=sqrt(sum((p1-p0).^2));
xy(2,1)=min(p1(1)*l/h,bx(2));
xy(2,2)=min(p1(2)*l/h,bx(4));
return
elseif npts == 3,
x=bx(2)-bx(1);
y=bx(4)-bx(3);
c=kd(1,2);
b=kd(1,3);
a=kd(2,3);
a2=a*a;
b2=b*b;
c2=c*c;
if max([a b c])<min(x,y),
if x>y,
if a>max(b,c),
xy(2,:)=[bx(1) bx(3)];
xy(3,:)=[bx(1)+a bx(3)];
cs=min(1,(a2+c2-b2)/2/a/c);
xy(1,1)=min(c*cs+bx(1),bx(2));
xy(1,2)=min(c*sqrt(1-cs*cs)+bx(3),bx(4));
return
elseif b>max(a,c),
xy(3,:)=[bx(1) bx(3)];
xy(1,:)=[bx(1)+b bx(3)];
cs=min(1,(b2+a2-c2)/2/b/a);
xy(2,1)=min(a*cs+bx(1),bx(2));
xy(2,2)=min(a*sqrt(1-cs*cs)+bx(3),bx(4));
return
else
xy(1,:)=[bx(1) bx(3)];
xy(2,:)=[bx(1)+c bx(3)];
cs=min(1,(b2+c2-a2)/2/b/c);
xy(3,1)=min(b*cs+bx(1),bx(2));
xy(3,2)=min(b*sqrt(1-cs*cs)+bx(3),bx(4));
return
end
else
if a>max(b,c),
xy(2,:)=[bx(1) bx(3)];
xy(3,:)=[bx(1) bx(3)+a];
cs=min((a2+c2-b2)/2/a/c,1);
xy(1,2)=min(bx(4),c*cs+bx(3));
xy(1,1)=min(c*sqrt(1-cs*cs)+bx(1),bx(2));
return
elseif b>max(a,c),
xy(3,:)=[bx(1) bx(3)];
xy(1,:)=[bx(1) bx(3)+b];
cs=min(1,(b2+a2-c2)/2/b/a);
xy(2,2)=min(bx(4),a*cs+bx(3));
xy(2,1)=min(a*sqrt(1-cs*cs)+bx(1),bx(2));
return
else
xy(1,:)=[bx(1) bx(3)];
xy(2,:)=[bx(1) bx(3)+c];
cs=min(1,(b2+c2-a2)/2/b/c);
xy(3,2)=min(bx(4),b*cs+bx(3));
xy(3,1)=min(b*sqrt(1-cs*cs)+bx(1),bx(2));
return
end
end
end
end
p0=bx([1 3]);
dela=ones(1000,1);
dela(500)=400;
dela=fft(dela);
if npts>78,
xy = zeros(npts,2);
xy(:,2) = (bx(3) + bx(4))/2;
xy(:,1) = linspace(bx(1), bx(2), npts)';
return
end
nones=ones(npts,1);
if any(kd(:)>0)
dx=bx(2)-bx(1);
dy=bx(4)-bx(3);
D=dx+dy;
d=D/20;
while max(2,ceil(dx/d))*max(2,ceil(dy/d))<200
d=d/2;
end
x=0:d:dx;
if x(end)<dx
x(end+1)=dx;
end
y=0:d:dy;
if y(end)<dy
y(end+1)=dy;
end
[x,y]=meshgrid(x,y);
x=x(:).';
y=y(:).';
k = length(x);
[dum,p]=sort(-max(kd));
kd=kd(p,p);
M=9;
N0=min(M,npts);
xy = zeros(npts,2);
xy(1:N0,:) = solverl(kd(1:N0,1:N0),dx,dy);
for N=M+1:npts
xy(N,:) = solver2(kd(1:N,1:N),N,x,y,xy(1:N,:),k,d,dx,dy);
end
for N=1:N0
xy(N,:) = solver2(kd,N,x,y,xy,k,d,dx,dy);
end
i0=0;
for N=1:N0
x1=xy(:,1);
y1=xy(:,2);
[XY1,XY2] = meshgrid(x1,y1);
dist = sqrt((XY1 - XY1').^2 + (XY2 - XY2').^2);

% Calculate strain matrix
strainMatrix = dist - kd;
strainMatrix(kd < 0) = 0;
[m,i]=max(sum(abs(strainMatrix)));
if i==i0, break; else, i0=i; end
xy(i,:) = solver2(kd,i,x,y,xy,k,d,dx,dy);
end
[v,w]=sort(p);
xy=xy(w,:);
xy=xy+p0(nones,:);
else
xy = p0(nones,:);
end

function xy = solverl(kd,dx,dy)
npts = length(kd);
[I,J,v] = find(triu(kd)>0);
k = 1000; % change k according to memory capacity.
x = rand(npts,k)*dx+i*rand(npts,k)*dy;
[tv,id] = min(sum(abs(abs(x(I,:)-x(J,:))-...
repmat(kd((J-1)*npts+I),1,k))));
xy=[real(x(:,id)) imag(x(:,id))];

function xy = solver2(kd,N,x,y,xy0,k,d,dx,dy)
[I,dm1,dm2] = find(kd(:,N)>0); nn=length(I);
nnones=ones(1,nn);
kones=ones(1,k);
[tv,id] = min(sum(abs(sqrt((xy0(I,kones)-x(nnones,:)).^2+...
(xy0(I,kones+1)-y(nnones,:)).^2)-kd(I,N+kones-1))));
x=max(0,x(id)-d):d/5:min(dx,x(id)+d);
y=max(0,y(id)-d):d/5:min(dy,y(id)+d);
[x,y]=meshgrid(x,y);
x=x(:).';
y=y(:).';
k = length(x);
nnones=ones(1,nn);
kones=ones(1,k);
[tv,id] = min(sum(abs(sqrt((xy0(I,kones)-x(nnones,:)).^2+...
(xy0(I,kones+1)-y(nnones,:)).^2)-kd(I,N+kones-1))));
xy=[x(id) y(id)];
```