# Diffing "Stuning02" and "Stuning03"

 Title: Stuning02 Stuning03 Author: Stijn Helsen Stijn Helsen Submitted: 2002-05-23 16:01:03 UTC 2002-05-23 16:01:36 UTC Status: Passed Passed Score: 7807.59 7807.56 Result: Average strain = 2031.390 Average strain = 2031.372 CPU Time: 167.161 168.361 Code: ```function xy = solver(kd,bx) npts0 = length(kd); dx=bx(2); dy=bx(4); D=dx+1i*dy; pm=D/2; xy=zeros(1,npts0)+pm; d=0.1*max(dx,dy); x=0:d:dx; y=0:d:dy; [xx,yy] = ndgrid(x,1i*y); XX=xx(:)+yy(:); dr=abs(D); kd(find(kd>dr))=dr; mx=max(kd-eye(npts0)); ismart=find(mx>=0); if isempty(ismart)|dx+dy<1e-3 xy=zeros(npts0,2); return; elseif length(ismart)==2 xy=zeros(npts0,2); f=kd(ismart(1),ismart(2))/dr; xy0=[0;D*f]; xy(ismart,:)=[real(xy0) imag(xy0)]; return end npts=length(ismart); kd=kd(ismart,ismart); S=(kd-eye(npts))>=0; nX=length(XX); xOnes=ones(nX,1); xZeros=zeros(nX,1); sCut = 1; olds = 1.0e30; [dum,p]=sort(-max(kd)); S=S(p,p); kd=kd(p,p); [v,w]=sort(p); if npts<50 xxx=11; xxx=12; else xxx=10; end; for iter=1:xxx 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 [s,M]=findstrain(S,xy,kd); 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); olds = 1.0e50; oldxy = xy; alpha = 1; for iter=1:40 [s,M]=findstrain(S,xy,kd); if( iter > 5 & abs(s-olds) / (s+1) < 0.0001 ) | s < sCut xy = oldxy; break; end if( s > olds ) xy = oldxy; break; end oldxy = xy; olds = s; for indx = 1:npts I = find(S(:,index)>0); dX=xy(I)-xy(indx); Dx=mean(dX.*(1-kd(I,index)./(abs(dX)+eps))); % Dx=mean(dX./(abs(dX)+eps) .* M(I,indx)); xy(indx) = xy(indx) + alpha*Dx; xy(indx) = min(dx,max(0,real(xy(indx))))+1i*min(dy,max(0,imag(xy(indx)))); end end alpha = 0.1; ngrad = 40; obj = 1.0e20*ones(ngrad,1); xybig = zeros(npts,ngrad); [objnew,strainMatrix]=findstrain(S,xy,kd); obj(1) = objnew; xybig(:,1) = xy; gradient = zeros(npts,1); for ij=2:ngrad [s,sM,dM]=findstrain(S,xy,kd); gradient(:) = sum(dM./(abs(dM)+eps).*sign(sM)); xy = xy - alpha * gradient; xy = min(max(real(xy),0),dx) + j * min(max(imag(xy),0),dy); [objnew,strainMatrix]=findstrain(S,xy,kd); obj(ij) = objnew; xybig(:,ij) = xy; if( obj(ij) > obj(ij-1) ) xy = xybig(:,ij-1); alpha = alpha / 3; end [bestobj,index] = min(obj); if( ij > 3 & index < 2 ) | bestobj < sCut break end end [bestobj,index] = min(obj); xy = xybig(:,index); xy0 = zeros(npts0,1); xy0(ismart)=xy; xy=[real(xy0) imag(xy0)]; function [s,strainMatrix,dX]=findstrain(S,xy,kd) sI=find(S); strainMatrix=S; dX=S; [X1,X2]=meshgrid(xy); dX(sI)=X1(sI)-X2(sI); strainMatrix(sI) = abs(dX(sI))-kd(sI); s = sum(abs(strainMatrix(sI)));```