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

Cyclops

by Nicholas Howe

Status: Passed
Results: 28758 (cyc: 6, node: 1979)
CPU Time: 1.742
Score: 28833.5
Submitted at: 2012-11-01 06:19:38 UTC
Scored at: 2012-11-01 17:24:38 UTC

Current Rank: 1620th (Highest: 11th )
Basis for: Cyclops (diff)

Comments
Please login or create a profile.
Code
function xy = cyclops(a, xy0, wts)

% Place nodes in radial pattern

xy = cyclopsR(a, xy0);
dxy = round(wts*(xy-xy0)./sum(wts));
xy = round(bsxfun(@minus,xy,dxy));
end

function xy = cyclopsR(a, xy0)

nn = length(a);
deg = sum(a);
la = logical(a);
xy = zeros(size(xy0));
left = true(1,nn);
fill = zeros(1,nn);
seq = zeros(1,nn);
ind = 0;
ap = zeros(1,nn);
ap1 = zeros(1,nn);
ap2 = zeros(1,nn);
nr = zeros(1,nn);
misc = 0.0001;

% highest degree node is root
rt = find(deg==max(deg),1);
left(rt) = false;

% add level 1 nodes radially
lvl = find(a(rt,:));
lvl = pickOrder(lvl,la(lvl,lvl));
ring = ceil(deg(rt)/8);
ringn = 8*ring;
spot = round(linspace(0,ringn-1,numel(lvl)));
ap(lvl) = spot/ringn;
r = ring./cos(mod(2*pi*ap(lvl)+pi/4,pi/2)-pi/4);
nr(lvl) = ring;
ap1(lvl) = (ap(lvl)+[ap(lvl(end))-1 ap(lvl(1:end-1))])/2+misc;
ap2(lvl) = (ap(lvl)+[ap(lvl(2:end)) ap(lvl(1))+1])/2-misc;
xy(lvl,:) = round([r.*cos(2*pi*ap(lvl));r.*sin(2*pi*ap(lvl))])';
left(lvl) = false;

%gplot(a(~left,~left),xy(~left,:),'o-');
%axis equal
%waitforbuttonpress;
        
% add remaining nodes level by level
lvlb = lvl;
lvl = find(left&any(a(lvl,:),1));
while any(lvl)
    % place children of previous level
    for par = lvlb
        ch = find(left&a(par,:));
        if isempty(ch), continue; end;
        ch = pickOrder(ch,la(ch,ch));
        nch = numel(ch);
        ring = nr(par)+1;
        while (floor(8*ring*ap2(par))-ceil(8*ring*ap1(par))+1 < nch)
            ring = ring+1;
        end;
        ringn = 8*ring;
        nr(ch) = ring;
        ap(ch) = round(linspace(ceil(8*ring*ap1(par)),floor(8*ring*ap2(par)),nch))/(8*ring);
        ap1(ch) = [ap1(par) (ap(ch(1:end-1))+ap(ch(2:end)))/2]+misc*(ap2(par)-ap1(par));
        ap2(ch) = [ap1(ch(2:end)) ap2(par)]-misc*(ap2(par)-ap1(par));
        xy(ch,:) = ringPts(ring,ap(ch));
        left(ch) = false;
        
        %clf
        %gplot(a(~left,~left),xy(~left,:),'o-');
        %hold on
        %gplot(a([par ch],[par ch]),xy([par ch],:),'ro-');
        %hold off
        %axis equal
        %waitforbuttonpress;
    end;

    % prepare for next level
    lvlb = lvl;
    lvl = find(left&any(a(lvl,:),1));
end;

%gplot(a,xy,'o-');
%axis equal

if any(left)
    rxy = cyclopsR(a(left,left), xy0(left,:));
    rxy(:,1) = rxy(:,1)-min(rxy(:,1))+max(xy(~left,:))+1;
    xy(left,:) = rxy;
end;
end

function xy = ringPts(ring,ap)
r = ring./cos(mod(2*pi*ap+pi/4,pi/2)-pi/4);
s = sin(2*pi*ap);
c = cos(2*pi*ap);
d = abs(mod(ap*8*ring+ring,2*ring)-ring);
%d = sqrt(r.^2-ring.^2);
xy = [ring*sign(c)' ring*sign(s)'];
sd = abs(s)>abs(c);
xy(sd,1) = d(sd).*sign(c(sd));
xy(~sd,2) = d(~sd).*sign(s(~sd));
end

function id = pickOrder(id0,la)
nid = numel(id0);
id = zeros(1,nid);
nc = sum(la);
nf = zeros(1,nid);
left = true(1,nid);
for i = 1:nid
    s = nf-nc;
    j = find(left&(s==max(s(left))),1);
    id(i) = id0(j);
    nf(la(j,:)) = nf(la(j,:))+1;
    nc(la(j,:)) = nc(la(j,:))-1;
    left(j) = false;
end;
end

function N = numKnots(XY,A)

% This function finds the number of knots in a network of points XY whose
% connectivity is given by A.

% Consider only unique lines (A is symmetric with ones on diagonal)
A = triu(A,1);
nLines = nnz(A);

% Find xy coordinates for first and second point on each line
[p1i,p2i] = find(A);
x1 = repmat(XY(p1i,1),1,nLines);
y1 = repmat(XY(p1i,2),1,nLines);
x2 = repmat(XY(p2i,1),1,nLines);
y2 = repmat(XY(p2i,2),1,nLines);
x3 = repmat(XY(p1i,1)',nLines,1);
y3 = repmat(XY(p1i,2)',nLines,1);
x4 = repmat(XY(p2i,1)',nLines,1);
y4 = repmat(XY(p2i,2)',nLines,1);

pick = tril(true(nLines),-1);
x1 = x1(pick);
y1 = y1(pick);
x2 = x2(pick);
y2 = y2(pick);
x3 = x3(pick);
y3 = y3(pick);
x4 = x4(pick);
y4 = y4(pick);

N = sum(areIntersecting(x1,y1,x2,y2,x3,y3,x4,y4));
end

function bool = areIntersecting(x1,y1,x2,y2,x3,y3,x4,y4)
% Determine if two line segments intersect.
%
% line1 = [x1, y1, x2, y2]
% line2 = [x3, y3, x4, y4]
%
% This function is robust to floating point precision issues assuming all
% input coordinates are integer-valued. If all coordinates have absolute
% values that do not exceed x, then the arithmetic operations used in this
% function could, as an intermediate step, produce a value with magnitude
% as large as 8*x^4. For double-precision floating point, this implies
% that the following must be satisfied:
%
% 8*x^4 < 2^53
%
% or equivalently, x < 5792.
%
% http://en.wikipedia.org/wiki/Line-line_intersection

bool = false(1,numel(x1));
gen = true(numel(x1),1);

% Test to see if the lines share exactly one endpoint. If so, the lines
% intersect if either of the free points lies on the line segment formed by
% the other two points.

s13 = find((x1==x3)&(y1==y3));
%s13 = find((x1==x3)&(y1==y3)&((x2~=x4)|(y2~=y4)));
bool(s13) = (isPointOnSegment(x2(s13), y2(s13), x1(s13), y1(s13), x4(s13), y4(s13))|isPointOnSegment(x4(s13), y4(s13), x1(s13), y1(s13), x2(s13), y2(s13)));
gen(s13) = false;
s14 = find((x1==x4)&(y1==y4));
%s14 = find((x1==x4)&(y1==y4)&((x2~=x3)|(y2~=y3)));
bool(s14) = (isPointOnSegment(x2(s14), y2(s14), x1(s14), y1(s14), x3(s14), y3(s14))|isPointOnSegment(x3(s14), y3(s14), x1(s14), y1(s14), x2(s14), y2(s14)));
gen(s14) = false;
s23 = find((x2==x3)&(y2==y3));
%s23 = find((x2==x3)&(y2==y3)&((x1~=x4)|(y1~=y4)));
bool(s23) = (isPointOnSegment(x1(s23), y1(s23), x2(s23), y2(s23), x4(s23), y4(s23))|isPointOnSegment(x4(s23), y4(s23), x1(s23), y1(s23), x2(s23), y2(s23)));
gen(s23) = false;
s24 = find((x2==x4)&(y2==y4));
%s24 = find((x2==x4)&(y2==y4)&((x1~=x3)|(y1~=y3)));
bool(s24) = (isPointOnSegment(x1(s24), y1(s24), x2(s24), y2(s24), x3(s24), y3(s24))|isPointOnSegment(x3(s24), y3(s24), x1(s24), y1(s24), x2(s24), y2(s24)));
gen(s24) = false;

% Next we check for parallel and coincident lines. Parallel lines
% obviously don't intersect. Coincident lines intersect if an endpoint from
% one of the lines lies on the other segment.

ss = find(gen&haveSameSlope(x1, y1, x2, y2, x3, y3, x4, y4));
gen(ss) = false;
si = ss(haveSameIntercept(x1(ss), y1(ss), x2(ss), y2(ss), x3(ss), y3(ss), x4(ss), y4(ss)));
bool(si) = ((isPointBetween(x3(si), x4(si), x1(si), 1) & isPointBetween(y3(si), y4(si), y1(si), 1)) | ...
               (isPointBetween(x3(si), x4(si), x2(si), 1) & isPointBetween(y3(si), y4(si), y2(si), 1)) | ...
               (isPointBetween(x1(si), x2(si), x3(si), 1) & isPointBetween(y1(si), y2(si), y3(si), 1)) | ...
               (isPointBetween(x1(si), x2(si), x4(si), 1) & isPointBetween(y1(si), y2(si), y4(si), 1)));

% If we get this far, we're in the general case of two well-defined
% non-parallel lines. The lines formed by the two segments must intersect
% somewhere. Find this intersection point and determine if it lies on the
% segments.

gen = find(gen);

Px_n = (x1(gen).*y2(gen) - y1(gen).*x2(gen)).*(x3(gen) - x4(gen)) - (x1(gen) - x2(gen)).*(x3(gen).*y4(gen) - y3(gen).*x4(gen));
Px_d = (x1(gen) - x2(gen)).*(y3(gen) - y4(gen)) - (y1(gen) - y2(gen)).*(x3(gen) - x4(gen));

Py_n = (x1(gen).*y2(gen) - y1(gen).*x2(gen)).*(y3(gen) - y4(gen)) - (y1(gen) - y2(gen)).*(x3(gen).*y4(gen) - y3(gen).*x4(gen));
Py_d = (x1(gen) - x2(gen)).*(y3(gen) - y4(gen)) - (y1(gen) - y2(gen)).*(x3(gen) - x4(gen));

bool(gen) = (isPointBetween(x1(gen), x2(gen), Px_n, Px_d) & ...
    isPointBetween(y1(gen), y2(gen), Py_n, Py_d) & ...
    isPointBetween(x3(gen), x4(gen), Px_n, Px_d) & ...
    isPointBetween(y3(gen), y4(gen), Py_n, Py_d));

end

function bool = haveSameIntercept(x1, y1, x2, y2, x3, y3, x4, y4)
bool = (x4 - x3).*(y1.*(x2 - x1) - x1.*(y2 - y1)) == ...
       (x2 - x1).*(y3.*(x4 - x3) - x3.*(y4 - y3));
end

function bool = haveSameSlope(x1, y1, x2, y2, x3, y3, x4, y4)
bool = (y2 - y1).*(x4 - x3) == (x2 - x1).*(y4 - y3);
end

function bool = isPointOnSegment(x1, y1, x2, y2, x3, y3)
% Determine if point (x1, y1) is on the segment formed by (x2, y2) and (x3, y3).
bool = haveSameSlope(x2, y2, x1, y1, x1, y1, x3, y3) & ...
       isPointBetween(x2, x3, x1, 1) & ...
       isPointBetween(y2, y3, y1, 1);
end

function bool = isPointBetween(pt1, pt2, num, den)
% Determine if either of the following is satisfied:
% pt1 <= num/den <= pt2
% pt1 >= num/den >= pt2

bool = ((pt1.*den <= num) & (num <= pt2.*den)) | ...
       ((pt1.*den >= num) & (num >= pt2.*den));
end