vanilla

function xyOut = solver(a, xyIn, m)
% Vicente Parot
    n = numel(xyIn)/2;
    [ii jj] = meshgrid(1:n,1:n);
    pi = [ii(:) jj(:)];
    pc = ~~a(pi(:,1)+n*(pi(:,2)-1));
    xyOut = xyIn + ceil(rand(size(xyIn))*3)-2;
    for it = 1:15
        s = gs(xyOut);
        s = ceil(abs(s)).*sign(s);
        xyOut = xyOut + 2*s;
    end
    s = gs(xyOut);
    s = ceil(abs(s)).*sign(s);
    xyOut = xyOut + s;
    r = rl(xyOut); 
    while r
        ism = true;
        while ism
            rn = randn(1,2);
            shift = ceil(abs(rn)).*sign(rn);
            [~,ism] = ismember(xyOut(r(1),:) + shift,xyOut,'rows');
        end
        xyt1 = xyOut;
        xyt2 = xyOut;
        xyt1(r(1),:) = xyt1(r(1),:) + shift;
        xyt2(r(2),:) = xyt2(r(2),:) + shift;
        if norm(gs(xyt1)) < norm(gs(xyt2))
            xyOut = xyt1;
        else
            xyOut = xyt2;
        end 
        r = rl(xyOut);
    end
    function s = gs(xyLocal)
        pxd = xyLocal(pi(:,1),1)-xyLocal(pi(:,2),1);
        pyd = xyLocal(pi(:,1),2)-xyLocal(pi(:,2),2);
        axf = -reshape(1./sqrt(pxd.^2+pyd.^2).*pxd.^2./(pxd.^2+pyd.^2).*sign(pxd),n,n);
        ayf = -reshape(1./sqrt(pxd.^2+pyd.^2).*pyd.^2./(pxd.^2+pyd.^2).*sign(pyd),n,n);
        cxf = reshape(pc.*sqrt(pxd.^2+pyd.^2).*pxd.^2./(pxd.^2+pyd.^2).*sign(pxd),n,n);
        cyf = reshape(pc.*sqrt(pxd.^2+pyd.^2).*pyd.^2./(pxd.^2+pyd.^2).*sign(pyd),n,n);
        axf = max(min(axf,2),-2);
        ayf = max(min(ayf,2),-2);
        axf = sum(axf,2);
        ayf = sum(ayf,2);
        cxf(isnan(cxf)) = 0;
        cyf(isnan(cyf)) = 0;
        cxf = sum(cxf,2);
        cyf = sum(cyf,2);
        s =     [axf./m' ayf./m']/50;
        s = s + [cxf./m' cyf./m']/10;
    end
    function p = rl(xy)
        p = numel(xy)~=numel(unique(xy,'rows'));
        if p
            for k = 1:n
                [~,p] = ismember(xy(k,:),xy((k+1):n,:),'rows');
                if p
                    p = [k p+k];
                    break
                end
            end
        end
    end
end