Winner Paulo Uribe (turbf1)

Finish 2002-05-23 00:00:00 UTC

Truncate2_wnb

by Paulo Uribe

Status: Passed
Results: Average strain = 2050.612
CPU Time: 144.587
Score: 7844.29
Submitted at: 2002-05-21 19:16:15 UTC
Scored at: 2002-05-21 19:44:41 UTC

Current Rank: 358th
Based on: Truncate2 (diff)
Basis for: Truncate2_wnb_wvec (diff)

Comments
Please login or create a profile.
Code
function xy = solver(kd,bx)
npts = length(kd); 
xy = zeros(npts,1); 
if any(kd(:)>0)
if npts<2,
return
elseif npts<10
 freemol=sum(sum(kd<0));
 if (freemol/(npts^2))>.25,
[I,J,v] = find(triu(kd)>0);
k = 1000; % change k according to memory capacity.
x = rand(npts,k)*bx(2); 
y = rand(npts,k)*bx(4);
[tv,id] = min(sum(abs(sqrt((x(I,:)-x(J,:)).^2+(y(I,:)-y(J,:)).^2)-...
	repmat(kd((J-1)*npts+I),1,k))));
xy=[x(:,id) y(:,id)];

msum=1e12;
mn=1;
for n = 1:300
x = repmat(xy(:,1),1,npts);
y = repmat(xy(:,2),1,npts);
delx = x - x';
dely = y - y';
r = sqrt(delx.^2 + dely.^2);
strain = (kd - r).*(kd ~= -1);
doubletotalstrain = sum(abs(strain(:)));
csum(n)=doubletotalstrain;
if csum(n)<msum,
    mxy=xy;
    msum=csum(n);
    mn=n;
end
if (n-mn)>20, break, end
if doubletotalstrain == 0, return, end
strain = strain/doubletotalstrain;
xy(:,1) = xy(:,1) - sum(delx.*strain)';
xy(:,2) = xy(:,2) - sum(dely.*strain)';
xy(:,1) = min(max(xy(:,1),0),bx(2));
xy(:,2) = min(max(xy(:,2),0),bx(4));
end
xy=mxy;
return
end
end
    nones=ones(npts,1); 
    dr=sqrt(bx(2)^2+bx(4)^2);
    kd(kd>dr)=dr;
    [dum,p]=sort(-max(kd));
    kd=kd(p,p);
    S=((kd-eye(npts))>=0); % Remember positive position
    sI=find(S); % Positive index
    strainMatrix = S; % Create strainMatrix;
    mx=-dum(1)*1.2;
    dx=min(bx(2),mx); 
    dy=min(bx(4),mx);
    D=dx+dy;
    d=D/16;
    while max(2,ceil(dx/d))*max(2,ceil(dy/d))<60
        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);
    xy0=x(:)+i*y(:);
    k = length(xy0);
    M=fix(8+npts/6);
    N0=min(M,npts);
    xy(1:N0) = solverl(S(1:N0,1:N0),kd(1:N0,1:N0),dx,dy,N0);
    for N=(M+1):npts
        xy(N) = solver2(S(1:N,N),kd(1:N,N),N,xy0,...
        xy(1:N),k,d,dx,dy,4);
    end
    for N=1:N0
        xy(N) = solver2(S(:,N),kd(:,N),N,[xy(N);xy0],xy,k+1,d,dx,dy,5);
    end
    [XY1,XY2] = meshgrid(xy);
    strainMatrix(sI) = abs(abs(XY1(sI)-XY2(sI))-kd(sI));
    sV = sum(strainMatrix);
    sN = find(sV>(sum(sV)/npts));
    for i=sN,
        xy(i) = solver2(S(:,i),kd(:,i),i,[xy(i);xy0],...
                   xy,k+1,d,dx,dy,5);
    end
    [v,w]=sort(p);
    xy=xy(w,:);
end
xy=[real(xy) imag(xy)];

function xy = solverl(S,kd,dx,dy,npts) 
[I,J] = find(triu(S,1)); 
if isempty(I),
    xy = zeros(npts,1);
    return;
end
y = kd(find(triu(S,1)));
k = 450+npts*20; % change k according to memory capacity.
x = zeros(npts,k);
x(:) = rand(npts*k,2)*[dx;i*dy];
[tv,id] = min(sum(abs(abs(x(I,:)-x(J,:))-y(:,ones(1,k)))));
xy=x(:,id);
xy = freeforce(S,kd,xy,dx,dy,npts,20);

function xy = freeforce(S,kd,xy,dx,dy,npts,n)
sI=find(S);
if isempty(sI),
    xy = zeros(npts,1);
end
fM=S;
d=kd(sI);
[X1,X2]=meshgrid(xy); 
dX=X1(sI)-X2(sI);
strain0=sum(abs(abs(dX)-d));
xy0=xy;
for iC=1:n,
   % Direction of force between every point and every other point
   % * magnitude of strain
   % * whether or not S(i,j) contains a don't-care for xy(i) and xy(j)
   
   fM(sI) = dX.*(kd(sI)./(abs(dX)+eps)-1);
   forces = sum(fM,2); 
   xy = xy0 - .1*forces;
   xy=max(zeros(npts,1),real(xy))+j*max(zeros(npts,1),imag(xy));
   xy=min(dx*ones(npts,1),real(xy))+j*min(dy*ones(npts,1),imag(xy));
   [X1,X2]=meshgrid(xy); 
   dX=X1(sI)-X2(sI);
   strain = sum(abs(abs(dX)-kd(sI)));
   if strain > strain0,
       xy = xy0;
       break;
   end
   xy0 = xy;
   strain0 = strain;
end;
        

function xy1 = solver2(S,kd,N,xy,xy0,k,d,dx,dy,n) 
I = find(S); 
if isempty(I),
  xy1 = xy0(N);
  return;
end 
z = kd(:); 
kones=ones(1,k); 
[XY1,XY2]=meshgrid(xy,xy0(I));
[tv0,id]=min(sum(abs(abs(XY1-XY2)-z(I,kones)))); 
xy1=xy(id);
for m=1:n
    d0=d;
    d=d0/4;
    [tv,xy1]=finegrid(I,xy1,xy0,d0,d,dx,dy,z);
    if tv0 - tv < 1e-6, break, end
    tv0 = tv;
end

function [tv,xy1]=finegrid(I,xy,xy0,d0,d,dx,dy,z);
x1=real(xy); 
y1=imag(xy);
x=max(0,x1-d0):d:min(dx,x1+d0); 
y=max(0,y1-d0):d:min(dy,y1+d0); 
[x,y]=meshgrid(x,y);
xy=x+i*y; 
xy=xy(:); 
k = length(xy); 
kones=ones(1,k); 
[XY1,XY2]=meshgrid(xy,xy0(I));
[tv,id]=min(sum(abs(abs(XY1-XY2)-z(I,kones)))); 
xy1=xy(id);