Finish 2012-11-07 16:00:00 UTC

Cheeeese

by Raphaël Candelier

Status: Passed
Results: 982 (cyc: 9, node: 1571)
CPU Time: 53.554
Score: 1071.79
Submitted at: 2012-11-07 15:19:48 UTC
Scored at: 2012-11-07 17:46:21 UTC

Current Rank: 1st (Highest: 1st )
Based on: Munster (diff)

Comments
Adi Pamungkas
27 Jan 2013
congratulation for Raphael Candeller
Tanmay
14 Apr 2013
congratulations
shehbaz Ali
03 Jul 2013
congratulation.
Raj Sodhi
04 Jul 2013
Congrats! Would you be kind enough to provide some sort of theory of operation on this code? I tried reading through it, but there was a lot of linear algebra stuff for which I could not understand the motivation.
Raphaël Candelier
14 Jul 2013
@Raj: sure, but can you be more explicit on what part(s) of the code you would like explanations ?
JEYASUTHA SANKARARAMAN
28 Oct 2013
good work raphael candelier....
Irfan Turk
16 Nov 2013
congratulation.
arun negi
31 Jan 2014
congrat
Nader Hadji Ghanbari
08 Mar 2014
Congrats
Alexis
17 Jun 2014
Congrats :)
Nay Linard
22 Oct 2014
What was the question for this code? By the way, I like it. Great work.
Please login or create a profile.
Code
function xyOut = solver(a, xyIn, wts)
% In code we trust

N = length(a);
tmp = pdist(xyIn);
W = max(tmp(:))*sum(wts);
nbin = 10;
xtra = 10;

diag(sum(a));
[V,~] = eig(ans-a, ans);

for ss=1:3
    [c{ss} c{ss+3}] = solver_1(a, ss+1, N, V+rand(size(V))/25);
end
c{7} = solver_2(xyIn, N, V(:,2:3),0);

[p1,p2] = find(triu(a,1));
[pickr,pickc] = find(tril(true(nnz(triu(a,1))),-1));
parfor ss=1:7
    NN(ss) = gradeIt(p1,p2,pickr,pickc,c{ss});
end

[~, best] = min(NN);
tmp = c{best};
xy = bsxfun(@minus, tmp, round(wts*(tmp-xyIn)./sum(wts)));

% --- Slipknot

m1 = min(xy,[], 1);
m2 = max(xy,[], 1);
xbin = unique(round(linspace(m1(1)-xtra, m2(1)+xtra, nbin)));
ybin = unique(round(linspace(m1(2)-xtra, m2(2)+xtra, nbin)));
[Y, X] = meshgrid(ybin, xbin);

for i = 1:N
    
    D = sqrt((xy(i,1)-X).^2 + (xy(i,2)-Y).^2);
    nei = find(a(i,:));
    I = p1~=i & p2~=i;
    v = funknots(xy, xy(i,1), xy(i,2), nei, p1(I), p2(I));
    
    if v
        
        K = funknots(xy, X, Y, nei, p1(I), p2(I));
        M = K+D*wts(i)/W;
        [m, mi] = min(M(:));
        if m<v
            x0 = X(mi);
            y0 = Y(mi);
            if ismember([x0 y0], xy, 'rows')
                for gs = 1:N
                    [X_, Y_] = meshgrid(x0+(-gs:gs), y0+(-gs:gs));
                    xy_ = [X_(:) Y_(:)];
                    xy_(ismember(xy_, xy, 'rows'),:) = NaN;
                    if any(~isnan(xy_(:,1)))
                        tmp = (x0-xy_(:,1)).^2+(y0-xy_(:,2)).^2;
                        tmp(isnan(tmp)) = Inf;
                        [~,mi] = min(tmp);
                        x0 = xy_(mi,1);
                        y0 = xy_(mi,2);
                        break;
                    end
                end
            end
            xy(i,:) = [x0 y0];
        end
    end
end

xyOut = xy;

end

function [xyOut xyOut1] = solver_1(a, ss, N, V)
randn('seed',ss);
xyOut=V(:,2:3);
sN=sqrt(N);
for ncut=1:2
    Mbeta=((1-ss/15)*eye(N)+(ss/15)*bsxfun(@rdivide,a+eye(N),sum(a,2)+1))^10;
    for n1=1:(117+ss)
        xyOut=Mbeta*xyOut;
        xc = bsxfun(@minus,xyOut,sum(xyOut)/N);
        [c1,D]=svd(xc');
        xyOut=xc*c1*diag(sN./(.1+D([1 4])))*c1';
    end
    [i,j]=find(a>0.002+ss*5.1e-4);
    sum((xyOut(i,:)-xyOut(j,:)).^2,2);
    k=find(ans>(3.15 + ss*0.102)*mean(ans));
    if isempty(k), break; end
    a(i(k)+N*(j(k)-1))=0.002+ss*5.1e-4;
end

xyOutH=xyOut;
single_idx = find(sum(a)==1);
[singles_link, ~] = find(a(:,single_idx)==1);
xyOut(single_idx,:)=xyOut(singles_link,:);

for zsingle=1:2
    xyOut0=sqrt(N)*detrend(xyOut,'constant')*diag(1./std(xyOut,1,1));
    k=5;
    [sxyOut0,idxequal]=sortrows(round(15*xyOut0));
    idxequal=idxequal(all(~diff(sxyOut0,1,1),2));
    xyOut=round(xyOut0);
    while size(unique(xyOut,'rows'),1)~=N
        k=k*(1.0+ss*0.055);
        xyOut=xyOut0*k;
        xyOut(idxequal,:)=(xyOut0(idxequal,:)+randn(numel(idxequal),2)/(9+ss))*k;
        xyOut=round(xyOut);
    end
    if zsingle==1
        xyOut1=xyOut;
        xyOut=xyOutH;
    end
end
end

function xyOut = solver_2(xyIn, N, V,xyOut)
mult    = norm(max(xyIn)-min(xyIn))/norm(max(V)-min(V));
while size(unique(xyOut,'rows'),1)<N
    xyOut       = round(V*mult);
    [~,i]       = unique(xyOut,'rows');
    m           = setdiff(1:N,i);
    xyOut(m,:)  = xyOut(m,:) + randi([-1 1],length(m),2);
    mult        = mult*1.5;
end
end

function ans = gradeIt(p1i,p2i,pickr,pickc,XYnew)
X = [XYnew(p1i,:) XYnew(p2i,:)];
x1 = X(pickr,:);
x3 = X(pickc,:);
sum(areIntersecting(x1(:,1),x1(:,2),x1(:,3),x1(:,4),x3(:,1),x3(:,2),x3(:,3),x3(:,4)));
end

function ans = areIntersecting(x1,y1,x2,y2,x3,y3,x4,y4)
a = x1.*y2 - y1.*x2;
b = x3.*y4 - y3.*x4;
Px_n = a.*(x3 - x4) - b.*(x1 - x2);
Py_n = a.*(y3 - y4) - b.*(y1 - y2);
Pxy_d = (x1 - x2).*(y3 - y4) - (y1 - y2).*(x3 - x4);
isPointBetween(x1, x2, Px_n, Pxy_d) & ...
    isPointBetween(y1, y2, Py_n, Pxy_d) & ...
    isPointBetween(x3, x4, Px_n, Pxy_d) & ...
    isPointBetween(y3, y4, Py_n, Pxy_d);
end

function ans = isPointBetween(pt1, pt2, num, den)
(pt1.*den - num).*(num - pt2.*den) >= 0;
end

function K = funknots(xy, X, Y, Lb, Lc, Ld)
% Knots for nuts

bx = permute(xy(Lb,1), [2 3 1]);
by = permute(xy(Lb,2), [2 3 1]);
cx = permute(xy(Lc,1), [2 3 4 1]);
cy = permute(xy(Lc,2), [2 3 4 1]);
dx = permute(xy(Ld,1), [2 3 4 1]);
dy = permute(xy(Ld,2), [2 3 4 1]);
ax_cx = bsxfun(@minus, X, cx);
ay_cy = bsxfun(@minus, Y, cy);
bx_ax = bsxfun(@minus, bx, X);
by_ay = bsxfun(@minus, by, Y);
dx_cx = bsxfun(@minus, dx, cx);
dy_cy = bsxfun(@minus, dy, cy);
nr = bsxfun(@minus, bsxfun(@times,ay_cy,dx_cx), bsxfun(@times,ax_cx,dy_cy));
ns = bsxfun(@minus, bsxfun(@times,ay_cy,bx_ax), bsxfun(@times,ax_cx,by_ay));
dn = bsxfun(@minus, bsxfun(@times,bx_ax,dy_cy), bsxfun(@times,by_ay,dx_cx));

r = bsxfun(@rdivide, nr, dn);
s = bsxfun(@rdivide, ns, dn);
K = sum(sum(r>=0 & r<1 & s>=0 & s<=1, 3), 4);
end

function D = pdist(X)
N = size(X,1);
D = NaN(N*(N-1)/2,1);
for i = 1:N-1
    D((i-1)*(N-i/2)+1:i*(N-(i+1)/2)) = (X(i,1)-X(i+1:end,1)).^2 + (X(i,1)-X(i+1:end,1)).^2;
end
D = sqrt(D);
end