# Diffing "GL6" and "Rana60"

 Title: GL6 Rana60 Author: Anon Rainer Lange Submitted: 2002-05-23 15:10:01 UTC 2002-05-23 15:26:06 UTC Status: Passed Passed Score: 7809.92 7809.94 Result: Average strain = 2032.593 Average strain = 2032.608 CPU Time: 169.873 169.672 Code: ```function xy = solver(kd,bx) npts = length(kd); if ~any(kd(:)>0), xy = zeros(npts,2); return end S=(kd-eye(npts))>=0; npts0=npts; nzS = find(sum(S)>0); S=S(nzS,nzS); sI=find(S); kd=kd(nzS,nzS); npts = length(nzS); nP=max(sum(S)); mx=max(kd(:))*(1+0.002*nP/npts); dx=min(bx(2),mx); dy=min(bx(4),mx); D=dx+1i*dy; pm=D/2; xy=zeros(1,npts)+pm; dr=abs(D); kd(find(kd>dr))=dr; minL = min(kd(find(triu(S,1)))); sCut = 0.1; d=max([max(dx,dy)/10, 6*sCut/nP, minL/2]); x=0:d:dx; y=0:d:dy; [xx,yy] = ndgrid(x,1i*y); XX=xx(:)+yy(:); nX=length(XX); xOnes=ones(nX,1); xZeros=zeros(nX,1); start=1; iones=ones(1,npts); sCut = 1; olds = 1.0e30; [dum,p]=sort(-max(kd)); S=S(p,p); sI=find(S); dX=S; kd=kd(p,p); [v,w]=sort(p); for iter=1:11 for index = 1:npts I = find(S(:,index)); if ~isempty(I) st = sum(abs(abs(XX(:,ones(size(I)))-xy(xOnes,I))-kd(xZeros+index,I)),2); [null,minloc] = min(st); xy(index) = XX(minloc); end end [X1,X2]=meshgrid(xy); dX(sI)=X1(sI)-X2(sI); s = sum(abs(abs(dX(sI))-kd(sI))); if((abs(s-olds)/(s+1)) < 0.002) | s < sCut break; end olds = s; end xy=xy.'; xy=xy(w,:); kd=kd(w,w); S=S(w,w); sI=find(S); dX=S; alpha = 0.1; ngrad = 60; obj = 1.0e20*ones(ngrad,1); xybig = zeros(npts,ngrad); gradient = zeros(npts,1); [X1,X2]=meshgrid(xy); dM=S;sM=S; dM(sI)=X1(sI)-X2(sI); sM(sI)=abs(dM(sI))-kd(sI); objnew = sum(abs(sM(sI))); gradient(:) = sum(dM./(abs(dM)+eps).*sign(sM)); ograd=gradient; obj(1) = objnew; xybig(:,1) = xy; for ij=2:ngrad xy = xy - alpha * gradient; xy = min(max(real(xy),0),dx) + j * min(max(imag(xy),0),dy); [X1,X2]=meshgrid(xy); dM=S;sM=S; dM(sI)=X1(sI)-X2(sI); sM(sI)=abs(dM(sI))-kd(sI); s = sum(abs(sM(sI))); gradient(:) = sum(dM./(abs(dM)+eps).*sign(sM)); obj(ij) = s; xybig(:,ij) = xy; if( obj(ij) > obj(ij-1) ) xy = xybig(:,ij-1); alpha = alpha / 10; gradient=ograd; else alpha = alpha * 2; alpha = alpha * 1.875; end [bestobj,index] = min(obj); if bestobj < sCut break end ograd = gradient; end [bestobj,index] = min(obj); xy = xybig(:,index); xy=[xy(:);zeros(npts0-npts,1)]; xy=[real(xy) imag(xy)];```