function [ X, m ] = A264041( n, varargin )
%A264041 Use cplex to get optimal configuration on n x n grid
%  Returns optimal configuration X and number of occupied cells m
%  If optional second argument is present, produces graphical output 
%  To produce the sequence A264041(1..N), use
%     A = zeros(1,N);
%     for n = 1:N
%       [~,A(n)] = A264041(n);
%     end
%     A
% 
nvars = 2*n^2;
% first n^2 constraints: x(ro,col,1) + x(ro,col,2) <= 1
Ai = [];
Aj = [];
conscount = 0;
for ro=1:n
    for col=1:n
        conscount = conscount+1;
        Ai(end+[1,2]) = conscount;
        Aj(end+[1,2]) = [varnum(ro,col,1,n),varnum(ro,col,2,n)];
    end
end
% next constraints: horizontal
for ro = 1:n
    for col=1:n-1
        conscount = conscount+2;
        Ai(end+[1:4]) = [conscount-1,conscount-1,conscount,conscount];
        Aj(end+[1:4]) = [varnum(ro,col,1,n),varnum(ro,col+1,2,n),varnum(ro,col,2,n),varnum(ro,col+1,1,n)];
    end
end
% next constraints: vertical
for ro = 1:n-1
    for col=1:n
        conscount = conscount+2;
        Ai(end+[1:4]) = [conscount-1,conscount-1,conscount,conscount];
        Aj(end+[1:4]) = [varnum(ro,col,1,n),varnum(ro+1,col,2,n),varnum(ro,col,2,n),varnum(ro+1,col,1,n)];
    end
end 
% next constraints: nw to se
for ro = 1:n-1
    for col = 1:n-1
       conscount = conscount+1;
       Ai(end+[1:2]) = conscount;
       Aj(end+[1:2]) = [varnum(ro,col,2,n),varnum(ro+1,col+1,2,n)];
    end
end
% last constraints: ne to sw
for ro = 1:n-1
    for col = 2:n
        conscount = conscount+1;
        Ai(end+[1:2]) = conscount;
        Aj(end+[1:2]) = [varnum(ro,col,1,n),varnum(ro+1,col-1,1,n)];
    end
end
As = ones(1,2*conscount);
Aineq = sparse(Ai,Aj,As);
bineq = ones(conscount,1);
obj = -ones(nvars,1);
[x,m, ~, output] = cplexbilp(obj,Aineq,bineq);
if output.cplexstatus ~= 101
    error(output.cplexstatusstring);
end
X = zeros(n,n);
m = -m;
for i=1:nvars
   if x(i) > 0.5
    [ro,col,p] = unpackvar(i,n);
    X(ro,col) = 1 - 2*(p-1);
   end
end
if nargin == 2
    fprintf('a(%d) = %d\n',n,m);
    disp(showconf(X));
end
end

function xn = varnum(ro,col,p,n)
% variable number for x(ro,col,p)
xn = (ro-1)*n + col + (p-1)*n^2;
end

function [ro,col,p] = unpackvar(xn,n)
% row, column, p for variable #n
if xn > n^2
    p = 2;
    xp = xn - n^2;
else
    p = 1;
    xp = xn;
end
col = mod(xp-1,n)+1;
ro = (xp - col)/n+1;
end

function S = showconf(X)
  n = size(X,1);
  S = char(' '*ones(n,n));
  S(X==1) = '/';
  S(X==-1) = '\';
end