Prolog GCHQ 2015 Christmas puzzle solution

% GCHQ 2015 Christmas puzzle solver inspired by techniques described at
% https://chriskohlhepp.wordpress.com/reasoning-systems/gchq-2015-crypto-puzzle-full-analysis/
% discussed at https://news.ycombinator.com/item?id=13349068
% original puzzle at https://www.gchq.gov.uk/news-article/christmas-card-cryptographic-twist-charity
%
% to run in SWI-Prolog: swipl -s gchq.pl -t 'preprocess(_)'
% output: a series of ever more refined QR codes ending with
%
% ███████ ███   █ █ ███████
% █     █ ██ ██     █     █
% █ ███ █     ███ █ █ ███ █
% █ ███ █ █  ██████ █ ███ █
% █ ███ █  █████ ██ █ ███ █
% █     █  ██       █     █
% ███████ █ █ █ █ █ ███████
%         ███   ███
% █ ██ ███  █ █ ███ █  █ ██
% █ █      ███ ██    █   █
%  ████ █ ████ ██ █    ██
%  █ █   █   █ █ ████ █ ███
%   ██  █ █ █      ██ █████
%    ███ ██ ██ ██████ ███ █
% █ █████████ █ █  ██    █
%  ██ █  ██   ██ ███     █
% ███ █ █ █  █    █████ █
%         █   ██ ██   █████
% ███████ █  ██   █ █ █ ███
% █     █ ██  █  ██   ██ █
% █ ███ █   ████  █████  █
% █ ███ █ ███ ██████████ ██
% █ ███ █ █  ██████ ██████
% █     █  ██      █ █ ██
% ███████ ██   █ ██   █████


% for transpose/2; on other Prologs than SWI-Prolog, implement
:- use_module(library(clpfd)).

row_constraints([[7,3,1,1,7],
                 [1,1,2,2,1,1],
                 [1,3,1,3,1,1,3,1],
                 [1,3,1,1,6,1,3,1],
                 [1,3,1,5,2,1,3,1],
                 [1,1,2,1,1],
                 [7,1,1,1,1,1,7],
                 [3,3],
                 [1,2,3,1,1,3,1,1,2],
                 [1,1,3,2,1,1],
                 [4,1,4,2,1,2],
                 [1,1,1,1,1,4,1,3],
                 [2,1,1,1,2,5],
                 [3,2,2,6,3,1],
                 [1,9,1,1,2,1],
                 [2,1,2,2,3,1],
                 [3,1,1,1,1,5,1],
                 [1,2,2,5],
                 [7,1,2,1,1,1,3],
                 [1,1,2,1,2,2,1],
                 [1,3,1,4,5,1],
                 [1,3,1,3,10,2],
                 [1,3,1,1,6,6],
                 [1,1,2,1,1,2],
                 [7,2,1,2,5]]).

col_constraints([[7,2,1,1,7],
                 [1,1,2,2,1,1],
                 [1,3,1,3,1,3,1,3,1],
                 [1,3,1,1,5,1,3,1],
                 [1,3,1,1,4,1,3,1],
                 [1,1,1,2,1,1],
                 [7,1,1,1,1,1,7],
                 [1,1,3],
                 [2,1,2,1,8,2,1],
                 [2,2,1,2,1,1,1,2],
                 [1,7,3,2,1],
                 [1,2,3,1,1,1,1,1],
                 [4,1,1,2,6],
                 [3,3,1,1,1,3,1],
                 [1,2,5,2,2],
                 [2,2,1,1,1,1,1,2,1],
                 [1,3,3,2,1,8,1],
                 [6,2,1],
                 [7,1,4,1,1,3],
                 [1,1,1,1,4],
                 [1,3,1,3,7,1],
                 [1,3,1,1,1,2,1,1,4],
                 [1,3,1,4,3,3],
                 [1,1,2,2,2,6,1],
                 [7,1,3,2,1,1]]).

input([[_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,█,█,_,_,_,_,_,_,_,█,█,_,_,_,_,_,_,_,█,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,█,█,_,_,█,_,_,_,█,█,_,_,█,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,█,_,_,_,_,█,_,_,_,_,█,_,_,_,█,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,█,█,_,_,_,_,█,█,_,_,_,_,█,_,_,_,_,█,█,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_],
       [_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_]]).

print_entry(Mark) :-
    (   atom(Mark)
    ->  write(Mark)
    ;   write('?') ).

print_line(Line) :-
    maplist(print_entry, Line),
    nl.

print_matrix(Matrix) :-
    maplist(print_line, Matrix).

% Generate a line of length TotalLength conforming to the given Constraints.
% Constraints is a list of numbers describing the lengths of blocks of marks
% on the line. For example, [3,2,1] means there is a block of three marks,
% one of two marks, and one of one mark. The rest of the line consists of
% blanks. Blocks of marks are separated by at least one blank, but there
% might be zero blanks before the first or after the last block.
% The total number of blanks on the line is tracked to limit unnecessary
% backtracking.
line(Constraints, TotalLength) -->
    { sum_list(Constraints, N) },
    { MaxInitialBlanks is TotalLength - N },
    { between(0, MaxInitialBlanks, InitialBlanks) },
    blanks(InitialBlanks),
    { BlankBudget is TotalLength - N - InitialBlanks },
    constraints(Constraints, BlankBudget, 0).

% Apply the constraints one by one. Each block of marks is followed by a
% block of at least one blank, except that the last block of marks may be
% followed by zero blanks.
constraints([C], Blanks, 0) -->
    marks(C),
    blanks(Blanks).
constraints([C,C1|Cs], BlankBudget, RemainingBlanks) -->
    marks(C),
    { between(1, BlankBudget, N) },
    blanks(N),
    { BlankBudget1 is BlankBudget - N },
    constraints([C1|Cs], BlankBudget1, RemainingBlanks).

blanks(0) -->
    [].
blanks(N) -->
    { N > 0 },
    [' '],
    { N1 is N - 1 },
    blanks(N1).

marks(0) -->
    [].
marks(N) -->
    { N > 0 },
    ['█'],
    { N1 is N - 1 },
    marks(N1).

% Apply the above DCG line//2.
constrain_line(Constraint, Line) :-
    phrase(line(Constraint, 25), Line).

constrain_line(Constraint, Line, Line) :-
    phrase(line(Constraint, 25), Line).

% Find all solutions of line//2 for the given Constraints and an initial
% line Givens which may contain previously placed marks.
constraint_lines(Constraint, Givens, Lines) :-
    setof(Line, Givens^constrain_line(Constraint, Givens, Line), Lines).

% Implement the interesting insight from the blog post linked above: All
% solutions to a given line constraint, and possibly starting from
% previously placed marks, may have something in common. For example, the
% N-th mark might be a mark in every solution, or the M-th mark might always
% be a blank. The indices of such always fixed "knowns" are collected in the
% Knowns structure, and immediately used to bind the corresponding elements
% of Givens.
knowns(Constraint, Givens, Knowns) :-
    constraint_lines(Constraint, Givens, Lines),
    transpose(Lines, Columns),
    known_marks_blanks(Columns, Knowns),
    bind_known_marks(Givens, Knowns).

known_marks_blanks(Columns, known_marks_blanks(Marks, Blanks)) :-
    known_marks_blanks(Columns, 1, [], [], MarksR, BlanksR),
    reverse(MarksR, Marks),
    reverse(BlanksR, Blanks).

known_marks_blanks([], _, Marks, Blanks, Marks, Blanks).
known_marks_blanks([Col|Cols], Idx, Marks0, Blanks0, Marks, Blanks) :-
    (   maplist(==('█'), Col)
    ->  Marks1 = [Idx|Marks0],
        Blanks1 = Blanks0
    ;   maplist(==(' '), Col)
    ->  Blanks1 = [Idx|Blanks0],
        Marks1 = Marks0
    ;   Marks1 = Marks0,
        Blanks1 = Blanks0 ),
    Idx1 is Idx + 1,
    known_marks_blanks(Cols, Idx1, Marks1, Blanks1, Marks, Blanks).

bind_known_marks(Line, known_marks_blanks(Marks, Blanks)) :-
    bind_known(Line, '█', 1, Marks),
    bind_known(Line, ' ', 1, Blanks).

bind_known([], _, _, _) :-
    !.
bind_known(_, _, _, []) :-
    !.
bind_known([E|Es], Mark, Idx, [M|Ms]) :-
    Idx1 is Idx + 1,
    (   M = Idx
    ->  E = Mark,
        bind_known(Es, Mark, Idx1, Ms)
    ;   bind_known(Es, Mark, Idx1, [M|Ms]) ).

% The preprocessing "derive-known-plaintext" operation from the blog post.
% As it turns out, this solves the whole puzzle without the need for a
% further backtracking search for solutions.
preprocess(Matrix) :-
    input(Matrix),
    % To experiment with all solutions to the given constraints, without
    % pre-marked fields (there are four such solutions), comment the
    % input(Matrix) line above and uncomment the following two lines that
    % set up an empty 25x25 matrix.
  % length(Matrix, 25),
  % maplist(same_length(Matrix), Matrix),
    transpose(Matrix, Transposed),
    row_constraints(RowConstraints),
    col_constraints(ColConstraints),
    preprocess(Matrix, Transposed, RowConstraints, ColConstraints, [], []).

% Iterate knowns/4 over the rows and columns of the matrix until a fixed
% point is reached.
preprocess(Matrix, Transposed, RowCs, ColCs, RowKs, ColKs) :-
    print_matrix(Matrix),
    nl,
    maplist(knowns, RowCs, Matrix, RowKs1),
    maplist(knowns, ColCs, Transposed, ColKs1),
    (   ( RowKs \= RowKs1 ; ColKs \= ColKs1 )
    ->  preprocess(Matrix, Transposed, RowCs, ColCs, RowKs1, ColKs1)
    ;   true ).