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);
|
